[Boost-commit] svn:boost r80596 - trunk/boost/polygon/detail

Subject: [Boost-commit] svn:boost r80596 - trunk/boost/polygon/detail
From: lucanus.j.simonson_at_[hidden]
Date: 2012-09-18 23:20:20

Author: ljsimons
Date: 2012-09-18 23:20:19 EDT (Tue, 18 Sep 2012)
New Revision: 80596
URL: http://svn.boost.org/trac/boost/changeset/80596

Log:
initial checkin of skeleton algorithm
Added:
   trunk/boost/polygon/detail/skeleton_detail.hpp   (contents, props changed)
Text files modified: 
   trunk/boost/polygon/detail/skeleton_predicates.hpp |     2 ++                                      
   trunk/boost/polygon/detail/voronoi_ctypes.hpp      |     9 +++++++++                               
   2 files changed, 11 insertions(+), 0 deletions(-)

Added: trunk/boost/polygon/detail/skeleton_detail.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/polygon/detail/skeleton_detail.hpp	2012-09-18 23:20:19 EDT (Tue, 18 Sep 2012)
@@ -0,0 +1,943 @@
+/*
+  Copyright 2012 Lucanus Simonson
+ 
+  Use, modification and distribution are subject to 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).
+*/
+#ifndef BOOST_POLYGON_SKELETON_DETAIL_HPP
+#define BOOST_POLYGON_SKELETON_DETAIL_HPP
+#include <boost/polygon/point_concept.hpp>
+#include <boost/polygon/segment_concept.hpp>
+#include <boost/polygon/polygon_set_concept.hpp>
+#include <boost/polygon/segment_utils.hpp>
+#include <queue>
+#include "skeleton_predicates.hpp"
+
+namespace boost { namespace polygon { 
+#ifdef BOOST_POLYGON_DEBUG1
+    inline void debug1(const char* msg) {std::cout << std::string(msg) << std::endl;}
+#else
+    inline void debug1(const char* msg) {}
+#endif
+
+
+    template <typename PointType, typename Segment2, typename Segment3>
+    bool
+    //computes in point the intersection of the lines of Segment2 and Segment3
+    //if Segment2 and Segment3 are parallel and non-colinear it returns false
+    //if segment2 and segment3 are parallel and colinear it returns the least low point in the direction
+    //both segments are directed, or if opposite directions it returns the midpoint between their two low points
+    generalized_intersection(PointType& point, const Segment2& segment2, const Segment3& segment3) {
+      typedef typename point_traits<PointType>::coordinate_type UnitOut;
+      UnitOut x11, x21, x12, x22, y12, y22, y11, y21, x_num, x_den, y_num, y_den, dx1, dy1, dx2, dy2;
+      x11 = UnitOut(segment2.low().get(HORIZONTAL));
+      x21 = UnitOut(segment3.low().get(HORIZONTAL));
+      y11 = UnitOut(segment2.low().get(VERTICAL));
+      y21 = UnitOut(segment3.low().get(VERTICAL));
+      x12 = UnitOut(segment2.high().get(HORIZONTAL));
+      x22 = UnitOut(segment3.high().get(HORIZONTAL));
+      y12 = UnitOut(segment2.high().get(VERTICAL));
+      y22 = UnitOut(segment3.high().get(VERTICAL));
+      dx1 = x12 - x11;
+      dy1 = y12 - y11;
+      dx2 = x22 - x21;
+      dy2 = y22 - y21;
+      if(dx1 * dy2 - dx2 * dy1 == 0) {
+        if((x22 - x11) * dx2 == (y22 - y11) * dy2)
+          return false; //lines have same slope but not colinear
+        UnitOut dll = (x11 - x21) * (x11 - x21) + (y11 - y21) * (y11 - y21);
+        UnitOut dlh = (x11 - x22) * (x11 - x22) + (y11 - y22) * (y11 - y22);
+        UnitOut dhl = (x12 - x21) * (x12 - x21) + (y12 - y21) * (y12 - y21);
+        UnitOut dhh = (x12 - x22) * (x12 - x22) + (y12 - y22) * (y12 - y22);
+        if(dhh < dhl && dhh < dlh && dhh < dll) {
+          //segments are directed toward eachother
+          point = PointType((x12+x22)/2,(y12+y22)/2); //average of high points
+        } else {
+          //away or in the same direction
+          point = PointType((x11+x21)/2,(y11+y21)/2); //average of low points
+        }
+      }
+      x_num = (x11 * dy1 * dx2 - x21 * dy2 * dx1 + y21 * dx1 * dx2 - y11 * dx1 * dx2);
+      x_den = (dy1 * dx2 - dy2 * dx1);
+      y_num = (y11 * dx1 * dy2 - y21 * dx2 * dy1 + x21 * dy1 * dy2 - x11 * dy1 * dy2);
+      y_den = (dx1 * dy2 - dx2 * dy1);
+      UnitOut x = x_num / x_den;
+      UnitOut y = y_num / y_den;
+      point = PointType(x, y);
+      std::cout << x11 << "," << y11 << "->" << x12 << "," << y12 << " " << x21 << "," << y21 << "->" << x22 << "," << y22 << "\n";
+      std::cout << "Intersection Point " << x << " " << y << std::endl;
+      return true;
+    }
+
+
+    template <typename Segment1, typename Segment2, typename Segment3>
+      bool
+    //computes in segment1 a segment with low point at the intersection of the lines
+    //of segment2 and segment3 oriented in the direction of the angle bisector of segment2
+    // and segment3
+    //degenerate case where segments are parallel and colinear
+    bisector(Segment1& segment1, const Segment2& segment2, const Segment3& segment3) {
+      typedef typename segment_traits<Segment1>::coordinate_type UnitOut;
+      point_data<UnitOut> gint;
+      if(generalized_intersection(gint, segment2, segment3)) {
+        UnitOut x11, x21, x12, x22, y12, y22, y11, y21, dx1, dy1, dx2, dy2;
+        segment1.low(gint);
+        x11 = UnitOut(segment2.low().get(HORIZONTAL));
+        x21 = UnitOut(segment3.low().get(HORIZONTAL));
+        y11 = UnitOut(segment2.low().get(VERTICAL));
+        y21 = UnitOut(segment3.low().get(VERTICAL));
+        x12 = UnitOut(segment2.high().get(HORIZONTAL));
+        x22 = UnitOut(segment3.high().get(HORIZONTAL));
+        y12 = UnitOut(segment2.high().get(VERTICAL));
+        y22 = UnitOut(segment3.high().get(VERTICAL));
+        dx1 = x12 - x11;
+        dy1 = y12 - y11;
+        dx2 = x22 - x21;
+        dy2 = y22 - y21;
+        
+        point_data<UnitOut> p1(x11, y11);
+        point_data<UnitOut> p2(x21, y21);
+        if(p1 == gint)
+          p1 = point_data<UnitOut>(x12, y12);
+        if(p2 == gint)
+          p2 = point_data<UnitOut>(x22, y22);
+        UnitOut dist1 = euclidean_distance(gint, p1);
+        UnitOut dist2 = euclidean_distance(gint, p2);
+        std::cout << dist1 << " " << dist2 << "\n";
+        UnitOut x_out = gint.x() + ((p1.x() - gint.x())*dist2/dist1 + (p2.x()-gint.x()))/2;
+        UnitOut y_out = gint.y() + ((p1.y() - gint.y())*dist2/dist1 + (p2.y()-gint.y()))/2;
+        point_data<UnitOut> end_point(x_out, y_out);
+        if(orientation(segment2, end_point) != orientation(segment3, end_point)) { 
+          std::cout << "we want the other bisector in this case\n";
+          y_out = gint.y() + ((p1.x() - gint.x())*dist2/dist1 + (p2.x()-gint.x()))/2;
+          x_out = gint.x() - ((p1.y() - gint.y())*dist2/dist1 + (p2.y()-gint.y()))/2;
+        }
+        end_point = point_data<UnitOut>(x_out, y_out);
+        segment1.high(end_point);
+        return true;
+      }
+      return false;
+    }
+
+    template <typename coordinate_type, typename output_coordinate_type>
+    struct skeleton_detail {
+
+      //DATA STRUCTURES
+      typedef segment_data<output_coordinate_type> edge_type;
+      typedef point_data<output_coordinate_type> vertex_type;
+      typedef segment_data<output_coordinate_type> segment_type;
+      
+      //forward declare
+      struct future_intersection;
+      struct skeleton_node;
+
+      //EVENT OBJECT can represent either split event or edge event with same type
+      struct event {
+        point_data<output_coordinate_type> point;
+        output_coordinate_type distance; //use for lazy compare, exact compare can be done with original coordinates in *parent
+        //event consists of intersection between angle bisectors of parent and parent->next or parent and split_on_next 
+        future_intersection* parent;
+        future_intersection* parent2;
+        //refer back to the edge between the target and its next
+        future_intersection* split_on_next; //if null is an edge event
+
+        bool operator<(const event& that) const {
+          debug1("event::operator<");
+          return distance < that.distance;
+        }
+        bool operator>(const event& that) const {
+          debug1("event::operator>");
+          return distance > that.distance;
+        }
+
+        output_coordinate_type x() const { return point.x(); }
+        void x(const output_coordinate_type& val) { point.x(val); }
+        
+        output_coordinate_type y() const { return point.y(); }
+        void y(const output_coordinate_type& val) { point.y(val); }
+        
+        output_coordinate_type r() const { return distance; }
+        void r(const output_coordinate_type& val) { distance = val; }
+      };
+      typedef typename boost::polygon::detail::skeleton_predicates
+      <typename boost::polygon::detail::skeleton_ctype_traits<coordinate_type> >::template event_formation_predicate<
+        segment_data<coordinate_type>, event> create_event;
+
+      //SET OF CIRCULAR LINKED LISTS OF FUTURE INTERSECTIONS
+      // each future intersection is associated with a pair of edges of the polygon
+      // originally coming from edges incident on each polygon vertex
+      // polygon holes are circular sub lists of each list such that split events can be found
+      struct future_intersection {
+        edge_type edge; //the edge that is next between this future intersection and the one stored in *next
+        future_intersection* prev;
+        future_intersection* next;
+        skeleton_node* parent_node;
+        output_coordinate_type active_distance;
+        //split events that refer to the edge between this and next
+        std::vector<typename std::set<event>::iterator> active_splits_on_next;
+        std::string label;
+        bool is_dead;
+        future_intersection() : prev(0x0), next(0x0), parent_node(0x0), is_dead(false), active_distance(-1.0) {
+          static char local_label = 'a';
+          label = local_label;
+          if(local_label == 'z')
+            local_label = 'a';
+          else
+            ++local_label;
+        }
+        ~future_intersection() {
+          debug1("~future_intersection()");
+          prev = 0x0;
+          next = 0x0;
+          parent_node = 0x0;
+        }
+      };
+
+      static void delete_future_intersection_ring(future_intersection* head) {
+        debug1("delete_future_intersection_ring()");
+        future_intersection* current = head;
+        if(head != 0x0)
+          if(head->prev != 0x0)
+            head->prev->next = 0x0; //give us a stopping point
+        while(current != 0x0) {
+          if(current->prev != 0x0)
+            current->prev->next = 0x0;
+          future_intersection* dead_node = current;
+          current = current->next;
+          delete dead_node;
+        }
+      }
+
+      //EVENT QUEUE of event objects by minimum weighted distance from edges
+      struct event_priority {
+        bool operator()(const event& e1, const event& e2) {
+          return e1 > e2;
+        }
+      };
+      typedef std::priority_queue<event, std::vector<event>, event_priority> event_queue_type;
+
+      //SKELETON OUTPUT DATA STRUCTURE
+      //OUTPUT REQUIREMENTS
+      //provide floating point representation of skeleton nodes and also references to all the input polygon edges from which 
+      //the node is weighted equidistant, that way original coordinates can be used by algorithms consuming skeleton data
+      //handle the case where there are more than three edges that are weighted cocirular by a zero length skeleton segment
+      //that can be collapsed in post processing
+      
+      struct skeleton_node {
+        skeleton_node(const point_data<output_coordinate_type>& pt) {
+          debug1("skeleton_node()");
+          static char local_label = 'A';
+          label = local_label;
+          if(local_label == 'Z')
+            local_label = 'A';
+          else
+            ++local_label;
+          vertex = pt;
+          linked_nodes[0] = linked_nodes[1] = linked_nodes[2] = 0x0;
+        }
+        void add_edge(skeleton_node* link) {
+          for(int i = 0; i < 3; ++i) {
+            if(!linked_nodes[i]) {
+              linked_nodes[i] = link;
+              break;
+            }
+          }
+        }
+        std::string label;
+        point_data<output_coordinate_type> vertex;
+        skeleton_node* linked_nodes[3]; //three is the common case
+        //for higher order nodes insert artificial zero length skeleton edges and duplicate node coordinates
+        //to create a 3-ary tree
+      };
+
+      //static member functions
+
+      //ORIGINAL POLYGON EDGES
+      //simply copy the original polygon edges into the future intersection edge pair data structure
+      static void initialize_skeleton_edges(std::vector<skeleton_node*>& skeleton, 
+                                            std::vector<future_intersection*>& all_faces,
+                                            event_queue_type& event_queue,
+                                            const std::vector<edge_type>& edges) {
+        debug1("initialize_skeleton_edges");
+        if(edges.empty()) return;
+        rectangle_data<output_coordinate_type> domain;
+        envelope_segments(domain, edges.begin(), edges.end());
+        future_intersection* head = 0x0, *current = 0x0, *prev = 0x0;
+        for(std::size_t i = 0; i < edges.size(); ++i) {
+          skeleton.push_back(new skeleton_node(point_data<output_coordinate_type>
+                                               (output_coordinate_type(edges[i].low().get(HORIZONTAL)),
+                                                output_coordinate_type(edges[i].low().get(VERTICAL)))));
+          current = new future_intersection();
+          all_faces.push_back(current);
+          if(i == 0) head = current;
+          current->edge = edges[i];
+          current->parent_node = (skeleton.back());
+          current->prev = prev;
+          current->is_dead = false;
+          if(prev) {
+            prev->next = current;
+          }
+          prev = current;
+        }
+        current->next = head;
+        head->prev = current;
+
+        current = head;
+        do {
+          debug1("loop");
+          if(orientation(current->prev->edge, current->edge) < 1) {
+            std::cout << "initialize reflex vertex\n";
+            future_intersection* inner = current->next->next;
+            while(inner != current) {
+              event e2;
+              if(compute_split_event(e2, current, inner)) {
+                if(e2.distance <= e2.parent->active_distance || e2.parent->active_distance < 0) {
+                  //split event has no parent2
+                  std::cout << "inserting split event\n";
+                  event_queue.push(e2);
+                  e2.parent->active_distance = e2.distance;
+                }
+              }
+              inner = inner->next;
+            }
+          }
+          event e;
+          if(compute_edge_event(e, current, current->next, domain)) {
+            if((e.distance <= e.parent->active_distance || e.parent->active_distance < 0) &&
+               (e.distance <= e.parent2->active_distance || e.parent2->active_distance < 0)) {
+              event_queue.push(e);
+              e.parent->active_distance = e.distance;
+              e.parent2->active_distance = e.distance;
+            }
+          }
+          current = current->next;
+        } while(current != head);
+      }
+
+      template<typename polygon_data_type>
+      static void initialize_skeleton_single_polygon(std::vector<skeleton_node*>& skeleton, 
+                                                     std::vector<future_intersection*>& all_faces,
+                                                     event_queue_type& event_queue,
+                                                     const polygon_data_type& polygon,
+                                                     bool reverse_winding = false) {
+        debug1("initialize_skeleton_single_polygon");
+        //polygons may be open or closed
+        if(polygon.size() == 0) return;
+        std::vector<edge_type> edges;
+        typename polygon_data_type::iterator_type itr = polygon.begin(), first, prev;
+        first = itr;
+        prev = itr;
+        for(++itr; itr != polygon.end(); ++itr) {
+          //ignore duplicate points
+          if(*prev != *itr) {
+            if(reverse_winding)
+              edges.push_back(edge_type(*itr, *prev));
+            else
+              edges.push_back(edge_type(*prev, *itr));
+            prev = itr;
+          }
+        }
+        if(*prev != *first){
+          if(reverse_winding)
+            edges.push_back(edge_type(*first, *prev));
+          else
+            edges.push_back(edge_type(*prev, *first));
+        }
+        if(reverse_winding) {
+          std::reverse(edges.begin(), edges.end());
+        }
+        initialize_skeleton_edges(skeleton, all_faces, event_queue, edges);
+      }
+
+      template<typename geometry_type>
+      static void initialize_skeleton(std::vector<skeleton_node*>& skeleton, std::vector<future_intersection*>& all_faces,
+                                      event_queue_type& event_queue,
+                                      const geometry_type& polygon_set) {
+        debug1("initialize_skeleton");
+        std::vector<polygon_with_holes_data<coordinate_type> > polys;
+        assign(polys, polygon_set); //this makes winding correct for polygons and holes
+        for(std::size_t i = 0; i < polys.size(); ++i) {
+          initialize_skeleton_single_polygon(skeleton, all_faces, event_queue, polys[i]);
+          for(typename polygon_with_holes_data<coordinate_type>::iterator_holes_type itrh = polys[i].begin_holes();
+              itrh != polys[i].end_holes(); ++itrh) {
+            initialize_skeleton_single_polygon(skeleton, all_faces, event_queue, *itrh);
+          }
+        }
+      }
+
+
+      //COMPUTE EDGE EVENT COORDINATES
+      //given four general line segments find the intersection between their angle bisectors or return false if it does not exist
+      //robustness concern: if two line segments don't share a vertex the bisector line will be anchored on their 
+      //projected intersection, which is a rational quantity, use the original coordinates of the line segments in the formula
+      //to compute the angle bisectors intersection rather than the anchor point of the bisector for computation
+      //perform quick predicates to check if the two bisectors should ever intersect
+      //include weight on each edge to bias the angle bisector
+      //include bound on intersection weighted distance 
+      //required degenerate case, if they are colinear then the bisector is perpendicular at their mid point between them on the line 
+      //(shared coordinate if abutting) if they are at 180 degrees
+      //optional degenerate case, if they are colinear but at zero degrees (spike) then the bisector is at the lesser of their 
+      //two start points and parallel to the segments
+
+      template <typename Segment>
+      static void print_segment(const Segment& s) {
+        std::cout << s.low().x() << "," << s.low().y() << "->" << s.high().x() << "," << s.high().y() << "\n";
+      }
+
+      static bool compute_split_event_point(point_data<output_coordinate_type>& point, 
+                                           edge_type e1, edge_type e2, edge_type e3, edge_type e4, edge_type e5) {
+        debug1("compute_split_event_point");
+        segment_type segment1, segment2, segment3, segment4, segment5;
+        print_segment(e1);
+        print_segment(e2);
+        print_segment(e4);
+        point_data<output_coordinate_type> v;
+        //compute that the intersection of the lines of e1 and e2 are on the correct side of the "opposite" segment e4
+        if(!generalized_intersection(v, e1, e2)) {
+          debug1("no intersection");
+          return false;
+        }
+        if(orientation(e4, v) <= 0) { //check that reflex vertex ray originates on the correct side of "opposite" edge e4
+          debug1("wrong side opposite");
+          return false;
+        }
+        //compute the bisector of e1 and e2, this is the reflex vertex bisector
+        if(bisector(segment1, e1, e2)) {
+          debug1("has bisector1");
+          if(bisector(segment2, e3, e4)) {
+            debug1("has bisector2");
+            if(bisector(segment3, e4, e5)) {
+              debug1("has bisector3");
+              //compute the edge event point using e1, e2, e2, e4 or e1, e2, e1, e4 if e4 is parallel to e2
+              if(bisector(segment4, e1, e4)) {
+                if(!generalized_intersection(point, segment1, segment4)) {
+                  debug1("no intersection 2");
+                  return false;
+                }
+              } else if(bisector(segment5, e2, e4)) {
+                if(!generalized_intersection(point, segment1, segment5)) {
+                  debug1("no intersection 2");
+                  return false;
+                }
+              } else {
+                debug1("no bisector4");
+                return false;
+              }
+              std::cout << point.x() << " " << point.y() << std::endl;
+              //test that point is on the coorect side of e1 and e2 and e4
+              if(orientation(e1, point) <= 0)
+                return false;
+              if(orientation(e2, point) <= 0)
+                return false;
+              if(orientation(e4, point) <= 0)
+                return false;
+              debug1("passed simple orientation tests");
+              //test that the point is on the correct side of the besector of e3,e4 and e4,e5
+              //it should be clockwise the bisector of e3,e4 and counterclockwise the bisector of e4,e5
+              //if e3,e4/e4e5 are concave then the correct sides will be reversed
+              if(orientation(e3, e4) == orientation(segment2, point))
+                return false;
+              debug1("passed complex orientation test1");
+              if(orientation(e4, e5) != orientation(segment3, point))
+                return false;
+              debug1("passed complex orientation test2");
+              //we seem to have a good split event point
+              return true;
+            }
+          }
+        }
+          
+        return false;
+      }
+
+      template <typename SegmentType>
+      static output_coordinate_type distance_squared_to_line(const point_data<output_coordinate_type>& pt, const SegmentType& segment) {
+        //(ax+by+c)/sqrt(a^2+b^2)
+        //(ax+by+c)(ax+by+c)/(a^2+b^2)
+        output_coordinate_type x11, x12, y12, y11, dx1, dy1, slope;
+        x11 = output_coordinate_type(segment.low().get(HORIZONTAL));
+        y11 = output_coordinate_type(segment.low().get(VERTICAL));
+        x12 = output_coordinate_type(segment.high().get(HORIZONTAL));
+        y12 = output_coordinate_type(segment.high().get(VERTICAL));
+        dx1 = x12 - x11;
+        dy1 = y12 - y11;
+        if(dx1 == 0.0)
+          return (pt.x() - x11)*(pt.x() - x11);
+        slope = dy1/dx1;
+        //0 = -slope * x + y + slope * x11 - y11
+        output_coordinate_type a = -slope;
+        // b = 1;
+        output_coordinate_type c = slope * x11 - y11;
+        output_coordinate_type numf = a*pt.x() + pt.y() + c;
+        return numf*numf/(a*a + 1);
+      }
+
+      static bool compute_split_event(event& e, future_intersection* first, future_intersection* opposite) {
+        debug1("compute_split_event");
+        point_data<output_coordinate_type> pt;
+        if(compute_split_event_point(e, first->prev->edge, first->edge, opposite->prev->edge, opposite->edge, opposite->next->edge)) {
+          std::cout << "found split event point\n";
+          segment_type s1;
+          assign(s1, opposite->edge);
+          print_segment(opposite->edge);
+          std::cout << pt.x() << " " << pt.y() << std::endl;
+          output_coordinate_type dist = distance_squared_to_line(pt, s1);
+          std::cout << "distance to split " << sqrt(dist) << std::endl;
+          segment_data<coordinate_type> segment1, segment2, segment3;
+          assign(segment1, first->prev->edge);
+          assign(segment2, first->edge);
+          assign(segment3, opposite->edge);
+          event e2;
+          create_event()(segment1, segment2, segment3, &e);
+          std::cout << "Andrii's point " << e.x() << " " << e.y() << " " << e.r() << std::endl;
+          e.point = pt;
+          e.distance = dist;
+          e.parent = first;
+          e.split_on_next = opposite;
+          return true;
+        }
+        return false;
+      }
+
+      static bool compute_edge_event_point(point_data<output_coordinate_type>& point, 
+                                           edge_type e1, edge_type e2, edge_type e3, edge_type e4) {
+        segment_type segment1, segment2;
+        if(bisector(segment1, e1, e2)) {
+          if(bisector(segment2, e3, e4)) {
+            return (generalized_intersection(point, segment1, segment2));
+          }
+        }
+        return false;
+          
+      }
+
+      static bool compute_edge_event(event& e, future_intersection* first, future_intersection* second, 
+                                     const rectangle_data<output_coordinate_type>& domain) {
+        point_data<output_coordinate_type> pt;
+        if(compute_edge_event_point(pt, first->prev->edge, first->edge, second->prev->edge, second->edge)) {
+          debug1("domain");
+          std::cout << pt.x() << " " << pt.y() << std::endl;
+          //edge event points outside (or on the boundary of) the polygon bounding box are useless and not worth processing
+          //they would never be reached by the forward progress of the algorithm through the event queue
+          //and would just slow it down and waste space
+          if(contains(domain, pt, false)) {
+            //compute distance from point to edge
+            segment_type s1;
+            assign(s1, first->edge);
+            output_coordinate_type dist = distance_squared_to_line(pt, s1);
+            output_coordinate_type dist2 = euclidean_distance(s1, pt);
+            std::cout << "distance " << dist2 << std::endl;
+            std::cout << "distance^2 " << dist << std::endl;
+            segment_data<coordinate_type> segment1, segment2, segment3;
+            assign(segment1, first->prev->edge);
+            assign(segment2, first->edge);
+            assign(segment3, second->edge);
+            event e2;
+            create_event()(segment1, segment2, segment3, &e);
+            std::cout << "Andrii's point " << e.x() << " " << e.y() << " " << e.r() << std::endl;
+            e.point = pt;
+            e.distance = dist;
+            e.parent = first;
+            e.parent2 = second;
+            e.split_on_next = 0x0;
+            return true;
+          }
+        }
+        return false;
+      }
+      
+      
+
+      //COMPUTE SPLIT EVENT COORDINATES
+      //given three general line segments find the intersection between any two angle bisectors of the three segments
+      //if two of the three segments are parallel choose the other angle bisector
+      //if all three are parallel (two should be colinear) find the point half way between the bisector of the colinear points
+      //the line of the third segment
+      
+      //TODO
+
+      static void find_edge_event(future_intersection* i1, future_intersection* i2,
+                                  event_queue_type& event_queue, const rectangle_data<output_coordinate_type>& domain) 
+      {
+        event new_event;
+        if(compute_edge_event(new_event, i1, i2, domain)) {
+          std::cout << "new event " << new_event.point.x() << " " << new_event.point.y() << std::endl;
+          if(orientation(i1->prev->edge, new_event.point) > 0 &&
+             orientation(i1->edge, new_event.point) > 0 &&
+             orientation(i2->edge, new_event.point) > 0) {
+            if((new_event.distance <= new_event.parent->active_distance || new_event.parent->active_distance < 0) &&
+               (new_event.distance <= new_event.parent->next->active_distance || new_event.parent->next->active_distance < 0)) {
+              event_queue.push(new_event);
+              new_event.parent->active_distance = new_event.distance;
+              new_event.parent->next->active_distance = new_event.distance;
+            } else {
+              std::cout << "distance test rejection " << new_event.parent->active_distance << "\n";
+            }
+          } else {
+            std::cout << "orientation test rejection\n";
+          }
+        }
+      }
+
+      static void print_ring(future_intersection* fi) {
+        future_intersection* current = fi;
+        do {
+          std::cout << current->label << "->";
+          current = current->next;
+        } while(current != fi);
+        std::cout << fi->label << std::endl;
+        do {
+          std::cout << current->label << "->";
+          current = current->prev;
+        } while(current != fi);
+        std::cout << fi->label << std::endl;
+      }
+
+      static void process_event(event current_event, std::vector<skeleton_node*>& skeleton, std::vector<future_intersection*>& all_faces,
+                                event_queue_type& event_queue, const rectangle_data<output_coordinate_type>& domain) {
+        if(current_event.split_on_next) {
+          //split_event
+          
+          std::cout << "split event " << current_event.point.x() << " " << current_event.point.y() << std::endl;
+          print_ring(current_event.parent);
+          //create one new skeleton node with the event parent future intersection's parent-node connected to it
+          skeleton.push_back(new skeleton_node(current_event.point));
+          current_event.parent->parent_node->add_edge(skeleton.back());
+          std::cout << current_event.parent->parent_node->label << "->" << skeleton.back()->label << std::endl;
+          skeleton.back()->add_edge(current_event.parent->parent_node);
+          current_event.parent->parent_node = skeleton.back();
+          //event parent node is not dead
+          //create one new future intersection associated with the with a copy of the "opposite" edge as its edge
+          //that split edge will appear in both rings
+          //the future intersection of the opposite edge gets linked to the parent of the split event as its new prev
+          //the newly created future intersection gets linked to the old prev of the parent of the split even as its new next
+          future_intersection* new_active = new future_intersection();
+          all_faces.push_back(new_active);
+          new_active->edge = current_event.split_on_next->edge;
+          new_active->prev = current_event.parent->prev;
+          new_active->next = current_event.split_on_next->next;
+          new_active->parent_node = skeleton.back();
+          new_active->active_distance = -1;
+          new_active->next->prev = new_active;
+          current_event.parent->prev->next = new_active;
+          current_event.parent->prev = current_event.split_on_next;
+          current_event.parent->prev = current_event.split_on_next;
+          current_event.parent->active_distance = -1;
+          current_event.split_on_next->next = current_event.parent;
+
+          find_edge_event(new_active, new_active->next, event_queue, domain);
+          find_edge_event(new_active->prev, new_active, event_queue, domain);
+
+          find_edge_event(current_event.split_on_next, current_event.split_on_next->next, event_queue, domain);
+          find_edge_event(current_event.split_on_next->next, current_event.split_on_next->next->next, event_queue, domain);
+          print_ring(current_event.split_on_next);
+          print_ring(new_active);
+        } else {
+          //edge event
+          std::cout << "event " << current_event.point.x() << " " << current_event.point.y() << std::endl;
+          std::cout << current_event.parent->parent_node->label << " " << current_event.parent2->parent_node->label << std::endl;
+          skeleton.push_back(new skeleton_node(current_event.point));
+          //this should implrement couter clockwise order
+          current_event.parent->parent_node->add_edge(skeleton.back());
+          std::cout << current_event.parent->parent_node->label << "->" << skeleton.back()->label << std::endl;
+          current_event.parent->next->parent_node->add_edge(skeleton.back());
+          std::cout << current_event.parent->next->parent_node->label << "->" << skeleton.back()->label << std::endl;
+          skeleton.back()->add_edge(current_event.parent->next->parent_node);
+          std::cout << skeleton.back()->label  << "->" << current_event.parent->next->parent_node->label << std::endl;
+          skeleton.back()->add_edge(current_event.parent->parent_node);
+          std::cout << skeleton.back()->label << "->" << current_event.parent->parent_node->label << std::endl;
+          //parent needs to be removed from the ring because it's face has collapsed at this event
+          current_event.parent->prev->next = current_event.parent->next;
+          current_event.parent->next->prev = current_event.parent->prev;
+          current_event.parent->is_dead = true;
+          current_event.parent->next->parent_node = skeleton.back();
+          if(current_event.parent->prev->next == current_event.parent->prev->prev) {
+            debug1("closed ring");
+            //we are done with this ring
+            current_event.parent->next->parent_node->add_edge(current_event.parent->prev->parent_node);
+            std::cout << current_event.parent->next->parent_node->label << "->" << current_event.parent->prev->parent_node->label << std::endl;
+            current_event.parent->prev->is_dead = true;
+            current_event.parent->prev->parent_node->linked_nodes[0] = current_event.parent->next->parent_node;
+            std::cout << current_event.parent->prev->parent_node->label << "->" << current_event.parent->next->parent_node->label << std::endl;
+            current_event.parent->next->is_dead = true;
+            return;
+          }
+          current_event.parent->active_distance = -1;
+          current_event.parent2->active_distance = -1;
+          //the following two blocks should be refactored
+          find_edge_event(current_event.parent->prev, current_event.parent->next, event_queue, domain);
+          find_edge_event(current_event.parent->next, current_event.parent->next->next, event_queue, domain);
+
+          //TODO, find split event here to handle ties between reflex vertices that generate a new reflect vertex
+        }
+      }
+
+      //main loop
+      static void execute_skeleton(std::vector<skeleton_node*>& skeleton, std::vector<future_intersection*>& all_faces,
+                                   event_queue_type& event_queue, const rectangle_data<output_coordinate_type>& domain, 
+                                   output_coordinate_type bound, bool handle_reflex_vertices = true) {
+        while(!event_queue.empty()) {
+          event current_event = event_queue.top();
+          event_queue.pop();
+          if(bound >=0 && current_event.distance > bound*bound)
+            break;
+          if(!current_event.parent->is_dead &&
+             (current_event.parent2 == 0x0 || !current_event.parent2->is_dead))
+            process_event(current_event, skeleton, all_faces, event_queue, domain);
+        }
+      }
+
+
+
+
+      //EAGERNESS
+      //it would be desirable when a split event occurs to recurse the algorithm into each sub-polygon to eagerly process it
+      //holes would need to be split between each sub-polygon and the event queue would need to be split as well
+      //the cost of splitting may more than offset the advantages
+      //if each future intersection of the sub polygon can be inserted into a new sub event queue they will be invalidated
+      //in the main event queue
+      //hole discovery could be costly 
+
+      //NOTE ABOUT RESIZING
+      //for skeleton based resizing all that is needed is to discard edges of the original polygon that collapse 
+      //with bounded edge events and handle reflex vertices of the orignal polygon that participate in bounded split events
+      //by associating their edges with the "opposite" split edge for computing bounded angle bisector output polygon vertex
+      //Simply terminate skeleton algorithm early and process the set of lists of active 
+
+      //BUILT IN TEST
+      template <typename stream_type>
+      static bool test_initialize_skeleton(stream_type& stdcout) {
+        debug1("test_initialize_skeleton");
+        std::string stdendl = "\n";
+        rectangle_data<coordinate_type> rect(100, 200, 300, 500);
+        rectangle_data<output_coordinate_type> domain(100, 200, 300, 500);
+        std::vector<skeleton_node*> skeleton; 
+        std::vector<future_intersection*> all_faces;
+        event_queue_type event_queue;
+        // initialize_skeleton(skeleton, all_faces, event_queue, rect);
+        // if(skeleton.size() == 4 && all_faces.size() == 4) {
+        //   stdcout << event_queue.size() << stdendl;
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) {
+        //     std::cout << (*itr) << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+        //       ((*itr)->linked_nodes[0]) << " " <<
+        //       ((*itr)->linked_nodes[1]) << " " << 
+        //       ((*itr)->linked_nodes[2]) << "\n";
+        //   }
+        //   execute_skeleton(skeleton, all_faces, event_queue, domain, 1000);
+        //   for(typename std::vector<future_intersection*>::iterator itr = all_faces.begin();
+        //       itr != all_faces.end(); ++itr)
+        //     delete *itr;
+        //   all_faces.clear();
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) {
+        //     std::cout << (*itr)->label << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+        //       ( ((*itr)->linked_nodes[0]) ? (*itr)->linked_nodes[0]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[1]) ? (*itr)->linked_nodes[1]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[2]) ? (*itr)->linked_nodes[2]->label : std::string("0"))  << "\n";
+        //   }
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) 
+        //     delete *itr;
+        //   skeleton.clear();
+        //   {
+        //   point_data<coordinate_type> pts[] = {point_data<coordinate_type>(120, 200),
+        //                                        point_data<coordinate_type>(300, 200),
+        //                                        point_data<coordinate_type>(300, 500),
+        //                                        point_data<coordinate_type>(100, 500),
+        //                                        point_data<coordinate_type>(100, 220),
+        //   };
+        //   polygon_data<coordinate_type> poly(pts, pts+5);
+        //   stdcout << "five point case\n";
+        //   initialize_skeleton(skeleton, all_faces, event_queue, poly);
+        //   execute_skeleton(skeleton, all_faces, event_queue, domain, 1000);
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) {
+        //     std::cout << (*itr)->label << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+        //       ( ((*itr)->linked_nodes[0]) ? (*itr)->linked_nodes[0]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[1]) ? (*itr)->linked_nodes[1]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[2]) ? (*itr)->linked_nodes[2]->label : std::string("0"))  << "\n";
+        //   }
+        //   }
+        //   for(typename std::vector<future_intersection*>::iterator itr = all_faces.begin();
+        //       itr != all_faces.end(); ++itr)
+        //     delete *itr;
+        //   all_faces.clear();
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) 
+        //     delete *itr;
+        //   skeleton.clear();
+
+        //   {
+        //   point_data<coordinate_type> pts[] = {point_data<coordinate_type>(120, 200),
+        //                                        point_data<coordinate_type>(280, 200),
+        //                                        point_data<coordinate_type>(300, 220),
+        //                                        point_data<coordinate_type>(300, 480),
+        //                                        point_data<coordinate_type>(280, 500),
+        //                                        point_data<coordinate_type>(120, 500),
+        //                                        point_data<coordinate_type>(100, 480),
+        //                                        point_data<coordinate_type>(100, 220),
+        //   };
+        //   polygon_data<coordinate_type> poly(pts, pts+8);
+        //   stdcout << "eight point case\n";
+        //   initialize_skeleton(skeleton, all_faces, event_queue, poly);
+        //   execute_skeleton(skeleton, all_faces, event_queue, domain, 1000);
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) {
+        //     std::cout << (*itr)->label << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+        //       ( ((*itr)->linked_nodes[0]) ? (*itr)->linked_nodes[0]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[1]) ? (*itr)->linked_nodes[1]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[2]) ? (*itr)->linked_nodes[2]->label : std::string("0"))  << "\n";
+        //   }
+        //   }
+
+        //   for(typename std::vector<future_intersection*>::iterator itr = all_faces.begin();
+        //       itr != all_faces.end(); ++itr)
+        //     delete *itr;
+        //   all_faces.clear();
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) 
+        //     delete *itr;
+        //   skeleton.clear();
+        //   {
+        //   point_data<coordinate_type> pts[] = {point_data<coordinate_type>(120, 200),
+        //                                        point_data<coordinate_type>(280, 200),
+        //                                        point_data<coordinate_type>(300, 220),
+        //                                        point_data<coordinate_type>(300, 480-100),
+        //                                        point_data<coordinate_type>(280, 500-100),
+        //                                        point_data<coordinate_type>(120, 500-100),
+        //                                        point_data<coordinate_type>(100, 480-100),
+        //                                        point_data<coordinate_type>(100, 220),
+        //   };
+        //   polygon_data<coordinate_type> poly(pts, pts+8);
+        //   stdcout << "co-circular eight point case\n";
+        //   initialize_skeleton(skeleton, all_faces, event_queue, poly);
+        //   execute_skeleton(skeleton, all_faces, event_queue, domain, 1000);
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) {
+        //     std::cout << (*itr)->label << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+        //       ( ((*itr)->linked_nodes[0]) ? (*itr)->linked_nodes[0]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[1]) ? (*itr)->linked_nodes[1]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[2]) ? (*itr)->linked_nodes[2]->label : std::string("0"))  << "\n";
+        //   }
+        //   }
+
+        //   for(typename std::vector<future_intersection*>::iterator itr = all_faces.begin();
+        //       itr != all_faces.end(); ++itr)
+        //     delete *itr;
+        //   all_faces.clear();
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) 
+        //     delete *itr;
+        //   skeleton.clear();
+
+
+
+        //   {
+        //   point_data<coordinate_type> pts[] = {point_data<coordinate_type>(100, 200),
+        //                                        point_data<coordinate_type>(150, 220),
+        //                                        point_data<coordinate_type>(200, 200),
+        //                                        point_data<coordinate_type>(180, 250),
+        //                                        point_data<coordinate_type>(200, 300),
+        //                                        point_data<coordinate_type>(150, 280),
+        //                                        point_data<coordinate_type>(100, 300),
+        //                                        point_data<coordinate_type>(120, 250),
+        //   };
+        //   polygon_data<coordinate_type> poly(pts, pts+8);
+        //   stdcout << "ninja star\n";
+        //   initialize_skeleton(skeleton, all_faces, event_queue, poly);
+        //   execute_skeleton(skeleton, all_faces, event_queue, domain, 1000);
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) {
+        //     std::cout << (*itr)->label << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+        //       ( ((*itr)->linked_nodes[0]) ? (*itr)->linked_nodes[0]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[1]) ? (*itr)->linked_nodes[1]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[2]) ? (*itr)->linked_nodes[2]->label : std::string("0"))  << "\n";
+        //   }
+        //   }
+
+        //   for(typename std::vector<future_intersection*>::iterator itr = all_faces.begin();
+        //       itr != all_faces.end(); ++itr)
+        //     delete *itr;
+        //   all_faces.clear();
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) 
+        //     delete *itr;
+        //   skeleton.clear();
+
+        //   {
+        //   point_data<coordinate_type> pts[] = {point_data<coordinate_type>(100, 200),
+        //                                        point_data<coordinate_type>(140, 240),
+        //                                        point_data<coordinate_type>(200, 200),
+        //                                        point_data<coordinate_type>(160, 240),
+        //                                        point_data<coordinate_type>(200, 300),
+        //                                        point_data<coordinate_type>(160, 260),
+        //                                        point_data<coordinate_type>(100, 300),
+        //                                        point_data<coordinate_type>(140, 260),
+        //   };
+        //   polygon_data<coordinate_type> poly(pts, pts+8);
+        //   stdcout << "evil ninja star\n";
+        //   initialize_skeleton(skeleton, all_faces, event_queue, poly);
+        //   std::cout << event_queue.size() << " GREP\n";
+        //   execute_skeleton(skeleton, all_faces, event_queue, domain, 1000);
+        //   for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+        //       itr != skeleton.end(); ++itr) {
+        //     std::cout << (*itr)->label << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+        //       ( ((*itr)->linked_nodes[0]) ? (*itr)->linked_nodes[0]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[1]) ? (*itr)->linked_nodes[1]->label : std::string("0"))  << " " <<
+        //       ( ((*itr)->linked_nodes[2]) ? (*itr)->linked_nodes[2]->label : std::string("0"))  << "\n";
+        //   }
+        //   }
+
+          for(typename std::vector<future_intersection*>::iterator itr = all_faces.begin();
+              itr != all_faces.end(); ++itr)
+            delete *itr;
+          all_faces.clear();
+          for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+              itr != skeleton.end(); ++itr) 
+            delete *itr;
+          skeleton.clear();
+
+          {
+          {
+          point_data<coordinate_type> pts[] = {point_data<coordinate_type>(100, 200),
+                                               point_data<coordinate_type>(300, 200),
+                                               point_data<coordinate_type>(300, 500),
+                                               point_data<coordinate_type>(100, 500),
+                                               point_data<coordinate_type>(290, 350),
+          };
+          polygon_data<coordinate_type> poly(pts, pts+5);
+          stdcout << "concave polygon\n";
+          initialize_skeleton(skeleton, all_faces, event_queue, poly);
+          execute_skeleton(skeleton, all_faces, event_queue, domain, 1000);
+          for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+              itr != skeleton.end(); ++itr) {
+            std::cout << (*itr)->label << " " << (*itr)->vertex.x() << " " << (*itr)->vertex.y() << " " <<
+              ( ((*itr)->linked_nodes[0]) ? (*itr)->linked_nodes[0]->label : std::string("0"))  << " " <<
+              ( ((*itr)->linked_nodes[1]) ? (*itr)->linked_nodes[1]->label : std::string("0"))  << " " <<
+              ( ((*itr)->linked_nodes[2]) ? (*itr)->linked_nodes[2]->label : std::string("0"))  << "\n";
+          }
+          }
+
+          for(typename std::vector<future_intersection*>::iterator itr = all_faces.begin();
+              itr != all_faces.end(); ++itr)
+            delete *itr;
+          all_faces.clear();
+          for(typename std::vector<skeleton_node*>::iterator itr = skeleton.begin();
+              itr != skeleton.end(); ++itr) 
+            delete *itr;
+          skeleton.clear();
+
+          return true;
+        }
+        return false;
+      }
+
+    };
+  }
+}
+
+#endif
+
+

Modified: trunk/boost/polygon/detail/skeleton_predicates.hpp
==============================================================================
--- trunk/boost/polygon/detail/skeleton_predicates.hpp	(original)
+++ trunk/boost/polygon/detail/skeleton_predicates.hpp	2012-09-18 23:20:19 EDT (Tue, 18 Sep 2012)
@@ -11,7 +11,9 @@
 #ifndef BOOST_POLYGON_DETAIL_SKELETON_PREDICATES
 #define BOOST_POLYGON_DETAIL_SKELETON_PREDICATES
 
+#ifndef BOOST_POLYGON_NO_DEPS
 #include <boost/cstdint.hpp>
+#endif
 
 #include <cmath>
 

Modified: trunk/boost/polygon/detail/voronoi_ctypes.hpp
==============================================================================
--- trunk/boost/polygon/detail/voronoi_ctypes.hpp	(original)
+++ trunk/boost/polygon/detail/voronoi_ctypes.hpp	2012-09-18 23:20:19 EDT (Tue, 18 Sep 2012)
@@ -10,7 +10,16 @@
 #ifndef BOOST_POLYGON_DETAIL_VORONOI_CTYPES
 #define BOOST_POLYGON_DETAIL_VORONOI_CTYPES
 
+#ifndef BOOST_POLYGON_NO_DEPS
 #include <boost/cstdint.hpp>
+#else
+namespace boost {
+  typedef int int32_t;
+  typedef long long int64_t;
+  typedef unsigned int uint32_t;
+  typedef unsigned long long uint64_t;
+};
+#endif
 
 #include <cmath>
 #include <cstring>