$include_dir="/home/hyper-archives/boost-commit/include"; include("$include_dir/msg-header.inc") ?>
Subject: [Boost-commit] svn:boost r78656 - in trunk/boost/geometry: algorithms/detail extensions/algorithms/buffer extensions/strategies
From: barend.gehrels_at_[hidden]
Date: 2012-05-26 16:58:58
Author: barendgehrels
Date: 2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
New Revision: 78656
URL: http://svn.boost.org/trac/boost/changeset/78656
Log:
[geometry] recent work on buffers (march/april) - pending commit
Text files modified: 
   trunk/boost/geometry/algorithms/detail/occupation_info.hpp                                  |   298 +++++++++++++++++-                      
   trunk/boost/geometry/extensions/algorithms/buffer/buffer_inserter.hpp                       |    21                                         
   trunk/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp                       |    38 ++                                      
   trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp             |   645 +++++++++++++++++++++++++++++++++++++-- 
   trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp |    44 +-                                      
   trunk/boost/geometry/extensions/algorithms/buffer/line_line_intersection.hpp                |     2                                         
   trunk/boost/geometry/extensions/algorithms/buffer/side_on_convex_range.hpp                  |     6                                         
   trunk/boost/geometry/extensions/strategies/buffer.hpp                                       |     9                                         
   8 files changed, 965 insertions(+), 98 deletions(-)
Modified: trunk/boost/geometry/algorithms/detail/occupation_info.hpp
==============================================================================
--- trunk/boost/geometry/algorithms/detail/occupation_info.hpp	(original)
+++ trunk/boost/geometry/algorithms/detail/occupation_info.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -9,6 +9,7 @@
 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OCCUPATION_INFO_HPP
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OCCUPATION_INFO_HPP
 
+ #define BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
 
 #include <algorithm>
 #include <boost/range.hpp>
@@ -18,6 +19,7 @@
 
 #include <boost/geometry/algorithms/equals.hpp>
 #include <boost/geometry/iterators/closing_iterator.hpp>
+#include <boost/geometry/iterators/ever_circling_iterator.hpp>
 
 
 namespace boost { namespace geometry
@@ -28,6 +30,53 @@
 namespace detail
 {
 
+template <typename P>
+class relaxed_less
+{
+    typedef typename geometry::coordinate_type<P>::type coordinate_type;
+
+    coordinate_type epsilon;
+
+public :
+
+    inline relaxed_less()
+    {
+        // TODO: adapt for ttmath, and maybe build the map in another way 
+        // (e.g. exact constellations of segment-id's, maybe adaptive.
+        //       for the moment take the sqrt (to consider it as comparable distances)
+        epsilon = std::numeric_limits<double>::epsilon() * 100.0;
+    }
+
+    inline bool operator()(P const& a, P const& b) const
+    {
+        coordinate_type const dx = math::abs(geometry::get<0>(a) - geometry::get<0>(b));
+        coordinate_type const dy = math::abs(geometry::get<1>(a) - geometry::get<1>(b));
+
+
+        if (dx < epsilon && dy < epsilon)
+        {
+            return false;
+        }
+        if (dx < epsilon)
+        {
+            return geometry::get<1>(a) < geometry::get<1>(b);
+        }
+
+        return geometry::get<0>(a) < geometry::get<0>(b);
+    }
+
+    inline bool equals(P const& a, P const& b) const
+    {
+        typedef typename geometry::coordinate_type<P>::type coordinate_type;
+
+        coordinate_type const dx = math::abs(geometry::get<0>(a) - geometry::get<0>(b));
+        coordinate_type const dy = math::abs(geometry::get<1>(a) - geometry::get<1>(b));
+
+        return dx < epsilon && dy < epsilon;
+    };
+};
+
+
 template <typename T, typename P1, typename P2>
 inline T calculate_angle(P1 const& from_point, P2 const& to_point)
 {
@@ -38,25 +87,35 @@
 }
 
 template <typename Iterator, typename Vector>
-inline Iterator advance_circular(Iterator it, Vector const& vector, bool forward = true)
+inline Iterator advance_circular(Iterator it, Vector const& vector, segment_identifier& seg_id, bool forward = true)
 {
         int const increment = forward ? 1 : -1;
         if (it == boost::begin(vector) && increment < 0)
         {
                 it = boost::end(vector);
+        seg_id.segment_index = boost::size(vector);
         }
         it += increment;
+    seg_id.segment_index += increment;
         if (it == boost::end(vector))
         {
+        seg_id.segment_index = 0;
                 it = boost::begin(vector);
         }
         return it;
 }
 
-template <typename T>
+template <typename Point, typename T>
 struct angle_info
 {
         typedef T angle_type;
+    typedef Point point_type;
+
+    segment_identifier seg_id;
+    int turn_index;
+    int operation_index;
+    Point intersection_point;
+    Point direction_point;
     T angle;
     bool incoming;
 };
@@ -79,7 +138,9 @@
                 }
         };
 
+public :
     collection_type angles;
+private :
     bool m_occupied;
         bool m_calculated;
 
@@ -115,12 +176,26 @@
                 , m_calculated(false)
     {}
 
-	template <typename Point1, typename Point2>
-	inline void add(Point1 const& point1, Point2 const& point2, bool incoming)
+	template <typename PointC, typename Point1, typename Point2>
+	inline void add(PointC const& map_point, Point1 const& direction_point, Point2 const& intersection_point,
+                    int turn_index, int operation_index,
+                    segment_identifier const& seg_id, bool incoming)
         {
+        //std::cout << "-> adding angle " << geometry::wkt(direction_point) << " .. " << geometry::wkt(intersection_point) << " " << int(incoming) << std::endl;
+		if (geometry::equals(direction_point, intersection_point))
+		{
+			//std::cout << "EQUAL! Skipping" << std::endl;
+			return;
+		}
+
         AngleInfo info;
         info.incoming = incoming;
-        info.angle = calculate_angle<typename AngleInfo::angle_type>(point1, point2);
+        info.angle = calculate_angle<typename AngleInfo::angle_type>(direction_point, map_point);
+        info.seg_id = seg_id;
+        info.turn_index = turn_index;
+        info.operation_index = operation_index;
+        info.intersection_point = intersection_point;
+        info.direction_point = direction_point;
         angles.push_back(info);
 
                 m_calculated = false;
@@ -136,12 +211,172 @@
                 return m_occupied;
         }
 
+    static inline void debug_left_turn(AngleInfo const& ai, AngleInfo const& previous)
+    {
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+		std::cout << "Angle: " << (ai.incoming ? "i" : "o") 
+            // << " " << si(ai.seg_id) 
+            << " " << (math::r2d * (ai.angle) )
+            << " turn: " << ai.turn_index << "[" << ai.operation_index << "]";
+
+        if (ai.turn_index != previous.turn_index
+            || ai.operation_index != previous.operation_index)
+        {
+            std::cout << " diff: " << math::r2d * math::abs(previous.angle - ai.angle);
+        }
+        std::cout << std::endl;
+#endif
+    }
+
+    static inline void debug_left_turn(std::string const& caption, AngleInfo const& ai, AngleInfo const& previous,
+                int code = -99, int code2 = -99, int code3 = -99, int code4 = -99)
+    {
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+		std::cout << " " << caption << ": " 
+                << (ai.incoming ? "i" : "o") 
+                << " " << (math::r2d * (ai.angle) )
+                << " turn: " << ai.turn_index << "[" << ai.operation_index << "]"
+                << " " << (previous.incoming ? "i" : "o") 
+                << " " << (math::r2d * (previous.angle) )
+                << " turn: " << previous.turn_index << "[" << previous.operation_index << "]";
+
+        if (code != -99)
+        {
+            std::cout << " code: " << code << " , " << code2 << " , " << code3 << " , " << code4;
+        }
+        std::cout << std::endl;
+#endif
+    }
+
+    template <typename Turns>
+    static inline void include_left_turn(AngleInfo const& outgoing, AngleInfo const& incoming,
+                    Turns const& turns,
+                    std::set<int> const& turn_indices, 
+                    std::set<std::pair<int, int> >& keep_indices)
+    {
+        // Safety check
+		if (turn_indices.count(outgoing.turn_index) > 0
+            && turn_indices.count(incoming.turn_index) > 0)
+		{
+            // Check if this is the intended turn (of the two segments making an open angle with each other)
+            if (turns[outgoing.turn_index].operations[outgoing.operation_index].other_id == 
+                turns[incoming.turn_index].operations[incoming.operation_index].seg_id)
+            {
+                //if (turns[outgoing.turn_index].operations[outgoing.operation_index].operation == detail::overlay::operation_union
+                //    || turns[outgoing.turn_index].operations[outgoing.operation_index].operation == detail::overlay::operation_continue)
+                {
+                    // This is the left turn. Mark it and don't check it for other pieces
+                    keep_indices.insert(std::make_pair(outgoing.turn_index, outgoing.operation_index));
+                    //std::cout << " INC";
+                }
+            }
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+            //else
+            //{
+            //    std::cout << " SKIP " << outgoing.turn_index << "[" << outgoing.operation_index << "]" << std::endl;
+            //}
+#endif
+        }
+    }
+
+
+    template <typename Turns>
+    inline void get_left_turns(Turns const& turns, 
+                    std::set<int> const& turn_indices, 
+                    std::set<std::pair<int, int> >& keep_indices)
+    {
+        typedef typename strategy::side::services::default_strategy
+		    <
+			    typename cs_tag<typename AngleInfo::point_type>::type
+		    >::type side_strategy;
+
+		std::sort(angles.begin(), angles.end(), angle_sort());
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+        std::cout << "get_left_turns()" << std::endl;
+#endif
+
+		std::size_t i = 0;
+		std::size_t n = boost::size(angles);
+
+		typedef geometry::ever_circling_range_iterator<collection_type const> circling_iterator;
+		circling_iterator cit(angles);
+
+        debug_left_turn(*cit, *cit);
+		for(circling_iterator prev = cit++; i < n; prev = cit++, i++)
+		{
+            debug_left_turn(*cit, *prev);
+
+            bool include = ! geometry::math::equals(prev->angle, cit->angle)
+				&& ! prev->incoming
+				&& cit->incoming;
+
+            if (include)
+			{
+                // Go back to possibly include more same left outgoing angles:
+                // Because we check on side too we can take a large "epsilon"
+                circling_iterator back = prev - 1;
+
+				typename AngleInfo::angle_type eps = 0.00001;
+                int b = 1;
+				for(std::size_t d = 0; 
+                        math::abs(prev->angle - back->angle) < eps 
+                            && ! back->incoming 
+                            && d < n;
+                    d++)
+				{
+                    --back;
+                    ++b;
+                }
+
+                // Same but forward to possibly include more incoming angles
+                int f = 1;
+                circling_iterator forward = cit + 1;
+				for(std::size_t d = 0;
+                        math::abs(cit->angle - forward->angle) < eps 
+                            && forward->incoming 
+                            && d < n;
+                        d++)
+				{
+					++forward;
+                    ++f;
+                }
+
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+                std::cout << "HANDLE " << b << "/" << f << " ANGLES" << std::endl;
+#endif
+
+		        for(circling_iterator bit = prev; bit != back; --bit)
+		        {
+                    int code = side_strategy::apply(bit->direction_point, prev->intersection_point, prev->direction_point);
+                    int code2 = side_strategy::apply(prev->direction_point, bit->intersection_point, bit->direction_point);
+		            for(circling_iterator fit = cit; fit != forward; ++fit)
+		            {
+                        int code3 = side_strategy::apply(fit->direction_point, cit->intersection_point, cit->direction_point);
+                        int code4 = side_strategy::apply(cit->direction_point, fit->intersection_point, fit->direction_point);
+
+                        if (! (code == -1 && code2 == 1))
+                        {
+                            include_left_turn(*bit, *fit, turns, turn_indices, keep_indices);
+                            debug_left_turn("ALSO", *fit, *bit, code, code2, code3, code4);
+                        }
+                    }
+                }
+			}
+		}
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+        std::cout << "USE " << keep_indices.size() << " OF " << turn_indices.size() << " TURNS" << std::endl;
+#endif
+    }
+
 };
 
 
 template <typename Point, typename Ring, typename Info>
-inline void add_incoming_and_outgoing_angles(Point const& point,
-				Ring const& ring, int segment_index,
+inline void add_incoming_and_outgoing_angles(Point const& map_point, Point const& intersection_point,
+				Ring const& ring, 
+                int turn_index,
+                int operation_index,
+                segment_identifier seg_id,
                 Info& info)
 {
     typedef typename boost::range_iterator
@@ -150,29 +385,42 @@
         >::type iterator_type;
 
         int const n = boost::size(ring);
-	if (segment_index >= n || segment_index < 0)
+	if (seg_id.segment_index >= n || seg_id.segment_index < 0)
         {
                 return;
         }
 
-	iterator_type it = boost::begin(ring) + segment_index;
+    segment_identifier real_seg_id = seg_id;
+	iterator_type it = boost::begin(ring) + seg_id.segment_index;
 
-    if (geometry::equals(point, *it))
+    // TODO: if we use turn-info ("to", "middle"), we know if to advance without resorting to equals
+    relaxed_less<Point> comparator;
+
+    if (comparator.equals(intersection_point, *it))
     {
-		it = advance_circular(it, ring, false);
+		// It should be equal only once. But otherwise we skip it (in "add")
+		it = advance_circular(it, ring, seg_id, false);
     }
 
-	info.add(*it, point, true);
+	info.add(map_point, *it, intersection_point, turn_index, operation_index, real_seg_id, true);
 
-	it = advance_circular(it, ring);
+    if (comparator.equals(intersection_point, *it))
+    {
+		it = advance_circular(it, ring, real_seg_id);
+	}
+	else
+	{
+		// Don't upgrade the ID
+		it = advance_circular(it, ring, seg_id);
+	}
     for (int defensive_check = 0; 
-		geometry::equals(point, *it) && defensive_check < n; 
+		comparator.equals(intersection_point, *it) && defensive_check < n; 
                 defensive_check++)
     {
-		it = advance_circular(it, ring);
+		it = advance_circular(it, ring, real_seg_id);
     }
 
-	info.add(*it, point, false);
+	info.add(map_point, *it, intersection_point, turn_index, operation_index, real_seg_id, false);
 }
 
 
@@ -183,10 +431,12 @@
 class occupation_map
 {
 public :
-    typedef std::map<Point, OccupationInfo, geometry::less<Point> > map_type;
+    typedef std::map<Point, OccupationInfo, relaxed_less<Point> > map_type;
+
     map_type map;
+	std::set<int> turn_indices;
 
-    inline OccupationInfo& find_or_insert(Point const& point)
+    inline OccupationInfo& find_or_insert(Point const& point, Point& mapped_point)
     {
         typename map_type::iterator it = map.find(point);
         if (it == boost::end(map))
@@ -195,7 +445,7 @@
                         = map.insert(std::make_pair(point, OccupationInfo()));
             it = pair.first;
         }
-
+        mapped_point = it->first;
         return it->second;
     }
 
@@ -204,6 +454,16 @@
         typename map_type::const_iterator it = map.find(point);
         return it != boost::end(map);
     }
+
+    inline bool contains_turn_index(int index) const
+    {
+        return turn_indices.count(index) > 0;
+    }
+
+    inline void insert_turn_index(int index)
+    {
+        turn_indices.insert(index);
+    }
 };
 
 
Modified: trunk/boost/geometry/extensions/algorithms/buffer/buffer_inserter.hpp
==============================================================================
--- trunk/boost/geometry/extensions/algorithms/buffer/buffer_inserter.hpp	(original)
+++ trunk/boost/geometry/extensions/algorithms/buffer/buffer_inserter.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -157,15 +157,16 @@
             output_point_type p;
             segment_type s1(previous_p1, previous_p2);
             segment_type s2(first_p1, first_p2);
-            line_line_intersection<output_point_type, segment_type>::apply(s1, s2, p);
-
-            std::vector<output_point_type> range_out;
-            join_strategy.apply(p, *begin, previous_p2, first_p1,
-                distance.apply(*(end - 1), *begin, side),
-                range_out);
-            if (! range_out.empty())
+            if (line_line_intersection<output_point_type, segment_type>::apply(s1, s2, p))
             {
-                collection.add_piece(buffered_join, *begin, range_out);
+                std::vector<output_point_type> range_out;
+                join_strategy.apply(p, *begin, previous_p2, first_p1,
+                    distance.apply(*(end - 1), *begin, side),
+                    range_out);
+                if (! range_out.empty())
+                {
+                    collection.add_piece(buffered_join, *begin, range_out);
+                }
             }
 
             // Buffer is closed automatically by last closing corner (NOT FOR OPEN POLYGONS - TODO)
@@ -367,8 +368,8 @@
 #ifdef BOOST_GEOMETRY_DEBUG_WITH_MAPPER
     //collection.map_offsetted(mapper);
         //collection.map_offsetted_points(mapper);
-    collection.map_turns(mapper);
-    //collection.map_opposite_locations(mapper);
+	//collection.map_turns(mapper);
+    collection.map_opposite_locations(mapper);
 #endif
 
     collection.discard_rings();
Modified: trunk/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp
==============================================================================
--- trunk/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp	(original)
+++ trunk/boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -55,7 +55,12 @@
                 state_type& state
                 )
     {
-std::cout << "WARNING " << reason << std::endl;
+#if defined(BOOST_GEOMETRY_COUNT_BACKTRACK_WARNINGS)
+extern int g_backtrack_warning_count;
+g_backtrack_warning_count++;
+#endif
+//std::cout << "!";
+//std::cout << "WARNING " << reason << std::endl;
 
         // TODO this is a copy of dissolve, check this for buffer
         state.m_good = false;
@@ -112,6 +117,11 @@
     
     int count_within, count_on_helper, count_on_offsetted, count_on_corner;
         int count_on_occupied;
+    int count_on_multi;
+#if defined(BOOST_GEOMETRY_COUNT_DOUBLE_UU)
+    int count_on_uu;
+#endif
+    std::set<int> piece_indices_to_skip;
     
 #ifdef BOOST_GEOMETRY_DEBUG_WITH_MAPPER
         std::string debug_string;
@@ -125,6 +135,10 @@
         , count_on_offsetted(0)
         , count_on_corner(0)
                 , count_on_occupied(0)
+        , count_on_multi(0)
+#if defined(BOOST_GEOMETRY_COUNT_DOUBLE_UU)
+        , count_on_uu(0)
+#endif
     {}
 };
 
@@ -133,6 +147,28 @@
 #endif // DOXYGEN_NO_DETAIL
 
 
+class si
+{
+private :
+    segment_identifier m_id;
+
+public :
+    inline si(segment_identifier const& id)
+        : m_id(id)
+    {}
+
+    template <typename Char, typename Traits>
+    inline friend std::basic_ostream<Char, Traits>& operator<<(
+            std::basic_ostream<Char, Traits>& os,
+            si const& s)
+    {
+        os << s.m_id.multi_index << "." << s.m_id.segment_index;
+        return os;
+    }
+};
+
+
+
 }} // namespace boost::geometry
 
 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_BUFFER_POLICIES_HPP
Modified: trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp
==============================================================================
--- trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp	(original)
+++ trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -20,8 +20,13 @@
 
 #include <boost/geometry/algorithms/equals.hpp>
 #include <boost/geometry/algorithms/covered_by.hpp>
+
 #include <boost/geometry/extensions/strategies/buffer_side.hpp>
 
+#include <boost/geometry/extensions/algorithms/buffer/buffered_ring.hpp>
+#include <boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp>
+#include <boost/geometry/extensions/algorithms/buffer/side_on_convex_range.hpp>
+
 #include <boost/geometry/algorithms/detail/overlay/calculate_distance_policy.hpp>
 #include <boost/geometry/algorithms/detail/overlay/enrichment_info.hpp>
 #include <boost/geometry/algorithms/detail/overlay/traversal_info.hpp>
@@ -31,11 +36,6 @@
 #include <boost/geometry/algorithms/detail/partition.hpp>
 
 
-#include <boost/geometry/extensions/algorithms/buffer/buffered_ring.hpp>
-#include <boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp>
-#include <boost/geometry/extensions/algorithms/buffer/side_on_convex_range.hpp>
-
-
 namespace boost { namespace geometry
 {
 
@@ -86,6 +86,60 @@
 };
 
 
+template <typename P>
+class relaxed_side
+{
+public :
+
+    // Template member function, because it is not always trivial
+    // or convenient to explicitly mention the typenames in the
+    // strategy-struct itself.
+
+    // Types can be all three different. Therefore it is
+    // not implemented (anymore) as "segment"
+
+    static inline int apply(P const& p1, P const& p2, P const& p)
+    {
+        typedef typename coordinate_type<P>::type coordinate_type;
+
+        coordinate_type const x = get<0>(p);
+        coordinate_type const y = get<1>(p);
+
+        coordinate_type const sx1 = get<0>(p1);
+        coordinate_type const sy1 = get<1>(p1);
+        coordinate_type const sx2 = get<0>(p2);
+        coordinate_type const sy2 = get<1>(p2);
+
+        // Promote float->double, small int->int
+        typedef typename geometry::select_most_precise
+            <
+                coordinate_type,
+                double
+            >::type promoted_type;
+
+        promoted_type const dx = sx2 - sx1;
+        promoted_type const dy = sy2 - sy1;
+        promoted_type const dpx = x - sx1;
+        promoted_type const dpy = y - sy1;
+
+        promoted_type const s 
+            = geometry::detail::determinant<promoted_type>
+                (
+                    dx, dy, 
+                    dpx, dpy
+                );
+
+        promoted_type const zero = promoted_type();
+        promoted_type const relaxed_epsilon = std::numeric_limits<double>::epsilon() * 5.0;
+
+        return math::abs(s) < relaxed_epsilon ? 0
+            : s > zero ? 1 
+            : -1;
+    }
+};
+
+
+
 template <typename Ring>
 struct buffered_piece_collection
 {
@@ -131,7 +185,7 @@
         // To check clustered locations we keep track of segments being opposite somewhere
         std::set<segment_identifier> m_in_opposite_segments;
 
-	struct buffer_occupation_info : public occupation_info<angle_info<coordinate_type> >
+	struct buffer_occupation_info : public occupation_info<angle_info<point_type, coordinate_type> >
         {
                 std::set<segment_identifier> seg_ids;
                 std::set<int> turn_indices;
@@ -243,11 +297,6 @@
 
                                 iterator next2 = next_point(ring2, it2);
 
-				if (piece1.index == 21 && piece2.index == 32)
-				{
-					int kkk = 0;
-				}
-
                 turn_policy::apply(*prev1, *it1, *next1,
                                     *prev2, *it2, *next2,
                                     the_model, std::back_inserter(m_turns));
@@ -294,20 +343,36 @@
         return segment_relation_disjoint;
         }
 
-    inline void add_angles(int index, point_type const& point, buffer_turn_operation<point_type> const& operation)
+    inline void add_angles(int turn_index, int operation_index, point_type const& point, buffer_turn_operation<point_type> const& operation)
     {
-        buffer_occupation_info& info = m_occupation_map.find_or_insert(point);
-        info.turn_indices.insert(index);
-        if (info.seg_ids.count(operation.seg_id) <= 0)
-        {
-            info.seg_ids.insert(operation.seg_id);
-            add_incoming_and_outgoing_angles(point, 
-						offsetted_rings[operation.seg_id.multi_index], 
-						operation.seg_id.segment_index, 
-						info);
-        }
+        point_type mapped_point;
+        buffer_occupation_info& info = m_occupation_map.find_or_insert(point, mapped_point);
+    	info.turn_indices.insert(turn_index);
+        info.seg_ids.insert(operation.seg_id);
+        add_incoming_and_outgoing_angles(mapped_point, point, 
+					offsetted_rings[operation.seg_id.multi_index], 
+                    turn_index, operation_index,
+					operation.seg_id, 
+					info);
     }
 
+    inline void add_angles(int turn_index)
+    {
+		if (! m_occupation_map.contains_turn_index(turn_index))
+		{
+			m_occupation_map.insert_turn_index(turn_index);
+
+			buffer_turn_info<point_type> const& turn = m_turns[turn_index];
+
+//std::cout << "Adding point " << turn_index << " " << geometry::wkt(turn.point) << std::endl;
+
+			add_angles(turn_index, 0, turn.point, turn.operations[0]);
+			add_angles(turn_index, 1, turn.point, turn.operations[1]);
+		}
+    }
+
+
+
     inline void classify_turn(buffer_turn_info<point_type>& turn, piece const& pc) const
     {
         if (pc.type == buffered_flat_end)
@@ -315,6 +380,12 @@
             return;
         }
 
+        // Don't check against a piece of which is was preferred in the "situations" for multi - map
+        if (turn.piece_indices_to_skip.count(pc.index) > 0)
+        {
+            return;
+        }
+
                 int flat_ends_involved = 0;
         for (int i = 0; i < boost::size(turn.operations); i++)
         {
@@ -346,11 +417,19 @@
                         return;
                 }
 
+        if (turn.operations[0].piece_index == 9 && turn.operations[1].piece_index == 35
+            && (pc.index == 3 || pc.index == 9 || pc.index == 35))
+        {
+            double epsilon = std::numeric_limits<double>::epsilon();
+            int kkk = 0;
+        }
+
+
                 segment_identifier on_segment_seg_id;
 
                 buffered_ring<Ring> const& ring = offsetted_rings[seg_id.multi_index];
 
-        int const side_offsetted = side_on_convex_range<side_strategy>(turn.point, 
+        int const side_offsetted = side_on_convex_range<side_strategy >(turn.point, 
                                                 boost::begin(ring) + seg_id.segment_index, 
                                                 boost::begin(ring) + pc.last_segment_index,
                                                 seg_id, on_segment_seg_id);
@@ -370,9 +449,6 @@
         }
         if (side_helper == 0)
         {
-			//BOOST_AUTO(d1, geometry::comparable_distance(turn.point, pc.helper_segments.back()));
-			//BOOST_AUTO(d2, geometry::comparable_distance(turn.point, pc.helper_segments.front()));
-            //if (d1 < 0.1 || d2 < 0.1)
                         if (geometry::equals(turn.point, pc.helper_segments.back())
                                 || geometry::equals(turn.point, pc.helper_segments.front()))
             {
@@ -397,43 +473,509 @@
         }
     }
 
+    inline void debug_preference(std::string const& caption, std::set<int> const& indices) const
+    {
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+        std::cout << caption << ": " << indices.size() << std::endl;
+        for (auto sit = indices.begin(); sit != indices.end(); ++sit)
+        {
+            int const index = *sit;
+            std::cout << "Keep "  << index // << "[" << sit->second << "]"
+                << " "<< si(m_turns[index].operations[0].seg_id)
+                << " "<< si(m_turns[index].operations[1].seg_id)
+                << " " << m_turns[index].operations[0].piece_index
+                << "/" << m_turns[index].operations[1].piece_index
+                << " " << method_char(m_turns[index].method)
+                << " " << operation_char(m_turns[index].operations[0].operation)
+                << "/" << operation_char(m_turns[index].operations[1].operation)
+				<< std::endl;
+        }
+#endif
+    }
+
+    template <typename T>
+    static inline std::pair<T, T> ordered_pair(T const& first, T const& second)
+    {
+        return first < second ? std::make_pair(first, second) : std::make_pair(second, first);
+    }
+
+    std::map<std::pair<segment_identifier, segment_identifier>, std::set<int> > turns_per_segment;
+    inline void fill_segment_map()
+    {
+        turns_per_segment.clear();
+        int index = 0;
+        for (typename boost::range_iterator<turn_vector_type>::type it =
+            boost::begin(m_turns); it != boost::end(m_turns); ++it, ++index)
+        {
+            turns_per_segment
+                [
+                    ordered_pair
+                        (
+                            m_turns[index].operations[0].seg_id,
+                            m_turns[index].operations[1].seg_id
+                        )
+                ].insert(index);
+        }
+    }
+
+    template <std::size_t Index, typename Turn>
+    static inline bool corresponds(Turn const& turn, segment_identifier const& seg_id)
+    {
+        return turn.operations[Index].operation == detail::overlay::operation_union
+            && turn.operations[Index].seg_id == seg_id;
+    }
+
+    inline bool prefer_one(std::set<int>& indices) const
+    {
+        std::map<segment_identifier, int> map;
+        for (auto sit = indices.begin(); sit != indices.end(); ++sit)
+        {
+            int const index = *sit;
+            map[m_turns[index].operations[0].seg_id]++;
+            map[m_turns[index].operations[1].seg_id]++;
+        }
+        std::set<segment_identifier> segment_occuring_once;
+        for (auto mit = map.begin(); mit != map.end(); ++mit)
+        {
+            if (mit->second == 1)
+            {
+                segment_occuring_once.insert(mit->first);
+            }
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_PREFER
+            std::cout << si(mit->first) << " " << mit->second << std::endl;
+#endif
+        }
+
+        if (segment_occuring_once.size() == 2)
+        {
+            // Try to find within all turns a turn with these two segments
+
+            std::set<segment_identifier>::const_iterator soo_it = segment_occuring_once.begin();
+            segment_identifier front = *soo_it;
+            soo_it++;
+            segment_identifier back = *soo_it;
+
+            std::pair<segment_identifier, segment_identifier> pair = ordered_pair(front, back);
+            auto it = turns_per_segment.find(pair);
+            if (it != turns_per_segment.end())
+            {
+                debug_preference("Found", it->second);
+                // Check which is the union/continue
+                segment_identifier good;
+                for (auto sit = it->second.begin(); sit != it->second.end(); ++sit)
+                {
+                    if (m_turns[*sit].operations[0].operation == detail::overlay::operation_union)
+                    {
+                        good = m_turns[*sit].operations[0].seg_id;
+                    }
+                    else if (m_turns[*sit].operations[1].operation == detail::overlay::operation_union)
+                    {
+                        good = m_turns[*sit].operations[1].seg_id;
+                    }
+                }
+
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_PREFER
+                std::cout << "Good: " << si(good) << std::endl;
+#endif
+
+                // Find in indexes-to-keep this segment with the union. Discard the other one
+                std::set<int> ok_indices;
+                for (auto sit = indices.begin(); sit != indices.end(); ++sit)
+                {
+                    if (corresponds<0>(m_turns[*sit], good) || corresponds<1>(m_turns[*sit], good))
+                    {
+                        ok_indices.insert(*sit);
+                    }
+                }
+                if (ok_indices.size() > 0 && ok_indices.size() < indices.size())
+                {
+                    indices = ok_indices;
+                    ////std::cout << "^";
+                    return true;
+                }
+            }
+        }
+        return false;
+    }
+
+	inline void get_left_turns(buffer_occupation_info& info)
+	{
+        std::set<std::pair<int, int> > keep_indices;
+		info.get_left_turns(m_turns, info.turn_indices, keep_indices);
+
+		if (keep_indices.empty())
+		{
+			return;
+		}
+
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+		for (std::set<std::pair<int, int> >::const_iterator sit = keep_indices.begin();
+			sit != keep_indices.end();
+			++sit)
+        {
+            std::cout << " " << sit->first << "[" << sit->second << "]";
+        }
+        std::cout << std::endl;
+#endif
+
+        // We have to create a version of the set with only turn_indices, ki
+        std::set<int> ki;
+		for (std::set<std::pair<int, int> >::const_iterator sit = keep_indices.begin();
+			sit != keep_indices.end();
+			++sit)
+		{
+            ki.insert(sit->first);
+
+            if (m_turns[sit->first].both(detail::overlay::operation_union))
+            {
+                // Avoid both turns of a u/u turn to be included.
+                if (keep_indices.count(std::make_pair(sit->first, 1 - sit->second)) == 0)
+                {
+                    m_turns[sit->first].operations[1 - sit->second].operation = detail::overlay::operation_none;
+                }
+#if defined(BOOST_GEOMETRY_COUNT_DOUBLE_UU)
+                else if (m_turns[sit->first].both(detail::overlay::operation_union))
+                {
+                    m_turns[sit->first].count_on_uu++;
+                }
+#endif
+            }
+
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+                std::cout << "Keep "  << sit->first << "[" << sit->second << "]"
+                    << " "<< si(m_turns[sit->first].operations[0].seg_id)
+                    << " "<< si(m_turns[sit->first].operations[1].seg_id)
+                    << " " << m_turns[sit->first].operations[0].piece_index
+                    << "/" << m_turns[sit->first].operations[1].piece_index
+                    << " " << method_char(m_turns[sit->first].method)
+                    << " " << operation_char(m_turns[sit->first].operations[0].operation)
+                    << "/" << operation_char(m_turns[sit->first].operations[1].operation)
+
+					<< std::endl;
+#endif
+
+        }
+
+        if (ki.size() == 2)
+        {
+            debug_preference("Keep", ki);
+            prefer_one(ki);
+            debug_preference("Kept", ki);
+        }
+
+		for (std::set<int>::const_iterator sit1 = info.turn_indices.begin();
+			sit1 != info.turn_indices.end();
+			++sit1)
+		{
+			if (ki.count(*sit1) == 0)
+			{
+        		m_turns[*sit1].count_on_multi++;
+            }
+            else
+            {
+		        for (std::set<int>::const_iterator sit2 = info.turn_indices.begin();
+			        sit2 != info.turn_indices.end();
+			        ++sit2)
+		        {
+			        if (sit2 != sit1)
+			        {
+				        m_turns[*sit1].piece_indices_to_skip.insert(m_turns[*sit2].operations[0].piece_index);
+				        m_turns[*sit1].piece_indices_to_skip.insert(m_turns[*sit2].operations[1].piece_index);
+			        }
+		        }
+            }
+        }
+	}
+
+	inline int piece_count(buffer_occupation_info const& info)
+	{
+		std::set<int> piece_indices;
+
+		for (std::set<int>::const_iterator sit = info.turn_indices.begin();
+			sit != info.turn_indices.end();
+			++sit)
+		{
+			piece_indices.insert(m_turns[*sit].operations[0].piece_index);
+			piece_indices.insert(m_turns[*sit].operations[1].piece_index);
+		}
+		return piece_indices.size();
+	}
+
         inline void classify_occupied_locations()
         {
+        fill_segment_map();
+
         for (typename boost::range_iterator<typename occupation_map_type::map_type>::type it =
                     boost::begin(m_occupation_map.map);
                         it != boost::end(m_occupation_map.map); ++it)
         {
                 buffer_occupation_info& info = it->second;
-			if (info.occupied())
+            int const pc = piece_count(info);
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+            std::cout << std::endl << "left turns: " << pc << " "
+                //<< std::setprecision(20)
+                << geometry::wkt(it->first) << std::endl;
+#endif
+			if (pc > 2)
                         {
-				for (std::set<int>::const_iterator sit = info.turn_indices.begin();
-					sit != info.turn_indices.end();
-					++sit)
+				if (info.occupied())
                                 {
-					m_turns[*sit].count_on_occupied++;
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION
+                    std::cout << "-> occupied" << std::endl;
+
+                    // std::set<int> turn_indices;
+                    //std::set<std::pair<int, int> > keep_indices;
+                    //info.get_left_turns(it->first, m_turns, turn_indices, keep_indices); // just for debug-info
+
+#endif
+					for (std::set<int>::const_iterator sit = info.turn_indices.begin();
+						sit != info.turn_indices.end();
+						++sit)
+					{
+						m_turns[*sit].count_on_occupied++;
+					}
                                 }
+				else
+				{
+					get_left_turns(info);
+				}
+	//std::cout << geometry::wkt(it->first) << " " << int(info.occupied()) << std::endl;
                         }
-//std::cout << geometry::wkt(it->first) << " " << int(info.occupied()) << std::endl;
         }
         }
 
+#define BOOST_GEOMETRY_DEBUG_BUFFER_SITUATION_MAP
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_SITUATION_MAP
+	inline int get_side(point_type const& point, Ring const& ring, int segment_index)
+	{
+		auto it = boost::begin(ring) + segment_index;
+		auto prev = it++;
+		return side_strategy::apply(point, *prev, *it);
+	}
+#endif
+
+    template <typename Iterator>
+    static inline point_type const& select_for_side(Iterator first, Iterator second, int index)
+    {
+        return index == 0 ? *first : *second;
+    }
+
+
+	inline int get_side(segment_identifier const& seg_id1, segment_identifier const& seg_id2, int which = 1) const
+	{
+        Ring const& ring1 = offsetted_rings[seg_id1.multi_index];
+        Ring const& ring2 = offsetted_rings[seg_id2.multi_index];
+
+		auto it1 = boost::begin(ring1) + seg_id1.segment_index;
+		auto it2 = boost::begin(ring2) + seg_id2.segment_index;
+
+		auto prev1 = it1++;
+		auto prev2 = it2++;
+
+		int code1 = side_strategy::apply(select_for_side(prev1, it1, which), *prev2, *it2);
+		int code2 = side_strategy::apply(select_for_side(prev2, it2, which), *prev1, *it1);
+
+        if (code1 == 1 && code2 == -1) return 1;
+        if (code1 == -1 && code2 == 1) return -1;
+
+        // ROBUSTNESS: in near collinear cases one might be zero, the other non-zero.
+        // This happens several times.
+        if (code1 != 0) return code1;
+        if (code2 != 0) return -code2;
+
+  //      // Check if the other side gives some more info
+  //      // (I've never seen this is the case though it might be so, if they are much longer.
+		//int code1f = side_strategy::apply(*prev1, *prev2, *it2);
+		//int code2f = side_strategy::apply(*prev2, *prev1, *it1);
+
+  //      if (code1f != 0 || code2f != 0)
+  //      {
+  //          std::cout << "From: " << code1f << " " << code2f << std::endl;
+  //          if (code1f != 0) return -code1f;
+  //          if (code2f != 0) return code2f;
+  //      }
+
+        // Collinear?
+#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_SITUATION_MAP
+        //std::cout << "Collinear: " << code1 << " " << code2 << std::endl;
+#endif
+        return 0;
+	}
+
+
+
+	inline void debug_segment(segment_identifier id)
+	{
+        auto const& ring = offsetted_rings[id.multi_index];
+		auto it = boost::begin(ring) + id.segment_index;
+		auto prev = it++;
+		geometry::model::referring_segment<point_type const&> segment(*prev, *it);
+        //std::cout << geometry::wkt(*prev) << " " << geometry::wkt(*it) << std::endl;
+	}
+
+
+	struct cluster_info
+	{
+		inline cluster_info(int i, point_type p, buffer_turn_operation<point_type> op)
+			: turn_index(i)
+			, point(p)
+			, operation(op)
+		{}
+
+		inline cluster_info()
+			: turn_index(-1)
+    	{}
+
+		int turn_index;
+		point_type point;
+		buffer_turn_operation<point_type> operation;
+	};
+
+    struct clustered_info
+    {
+        std::vector<cluster_info> intersecting_segments;
+        std::set<segment_identifier> intersecting_ids;
+        int piece_index;
+    };
+
+#ifdef OLD
+    struct situation_info
+    {
+        std::set<int> turn_indices;
+        std::set<segment_identifier> seg_ids;
+    };
+#endif
+
+	inline bool add_mutual_intersection(clustered_info const& cluster, segment_identifier const& seg_id)
+	{
+		bool result = false;
+        //if (cluster.intersecting_ids.count(seg_id) > 0)
+		for(auto it = cluster.intersecting_segments.begin(); it != cluster.intersecting_segments.end(); it++)
+		{
+			if (it->operation.seg_id == seg_id)
+			{
+                //add_angles(it->turn_index);
+				//add_angles(it->turn_index, it->point, it->operation);
+				result = true;
+    		}
+		}
+		return result;
+	}
+
+	inline void add_mutual_intersections(clustered_info const& cluster, std::map<segment_identifier, clustered_info> const& map)
+	{
+		for(auto it1 = cluster.intersecting_segments.begin(); it1 != cluster.intersecting_segments.end(); it1++)
+		{
+			auto const& other_cluster_it = map.find(it1->operation.seg_id);
+			if (other_cluster_it != map.end())
+			{
+				for(auto it2 = it1 + 1; it2 != cluster.intersecting_segments.end(); it2++)
+				{
+					if (! m_occupation_map.contains_turn_index(it1->turn_index)
+						|| ! m_occupation_map.contains_turn_index(it2->turn_index))
+					{
+						// Either these two segments intersect, or they are perfectly collinear.
+						// First check the intersection:
+						if (add_mutual_intersection(other_cluster_it->second, it2->operation.seg_id))
+						{
+							////std::cout << "#";
+							add_angles(it1->turn_index);
+							add_angles(it2->turn_index);
+						}
+						else
+						{
+							// Then the collinearity
+							int const code = get_side(it1->operation.seg_id, it2->operation.seg_id);
+							if (code == 0)
+							{
+								////std::cout << "1";
+								add_angles(it1->turn_index);
+								add_angles(it2->turn_index);
+							}
+							else 
+							{
+								int const code = get_side(it1->operation.seg_id, it2->operation.seg_id, 0);
+								if (code == 0)
+								{
+									////std::cout << "0";
+									add_angles(it1->turn_index);
+									add_angles(it2->turn_index);
+								}
+								else if (geometry::equals(it1->point, it2->point))
+								{
+									////std::cout << "*";
+									add_angles(it1->turn_index);
+									add_angles(it2->turn_index);
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+
+    inline void add_multi_intersections_to_occupation_map()
+	{
+        // Pass 1: create map of all segments
+		typedef std::map<segment_identifier, clustered_info> map_type;
+		map_type map;
+
+		int index = 0;
+        for (typename boost::range_iterator<turn_vector_type>::type it =
+            boost::begin(m_turns); it != boost::end(m_turns); ++it, ++index)
+        {
+			buffer_turn_info<point_type> const& turn = *it;
+
+			// if (! turn.both(detail::overlay::operation_union))
+			{
+
+				// Take care with all the indices
+				map[turn.operations[0].seg_id].piece_index = turn.operations[0].piece_index;
+				map[turn.operations[0].seg_id].intersecting_segments.push_back(cluster_info(index, turn.point, turn.operations[1]));
+				map[turn.operations[0].seg_id].intersecting_ids.insert(turn.operations[1].seg_id);
+
+				map[turn.operations[1].seg_id].piece_index = turn.operations[1].piece_index;
+				map[turn.operations[1].seg_id].intersecting_segments.push_back(cluster_info(index, turn.point, turn.operations[0]));
+				map[turn.operations[1].seg_id].intersecting_ids.insert(turn.operations[0].seg_id);
+			}
+		}
+
+        // Pass 2: 
+        // Verify all segments crossing with more than one segment, and if they intersect each other,
+        // at that pair
+		for (typename map_type::const_iterator mit = map.begin(); mit != map.end(); ++mit)
+		{
+			if (mit->second.intersecting_segments.size() > 1)
+			{
+				add_mutual_intersections(mit->second, map);
+			}
+		}
+    }
+
         inline void get_occupation()
     {
                 fill_opposite_segments();
 
+        // Pass 1: fill all segments part of opposite segments
         int index = 0;
         for (typename boost::range_iterator<turn_vector_type>::type it =
             boost::begin(m_turns); it != boost::end(m_turns); ++it, ++index)
         {
                         buffer_turn_info<point_type>& turn = *it;
+//std::cout << "Referring to point " << geometry::wkt(turn.point) << std::endl;
                         if (m_in_opposite_segments.count(turn.operations[0].seg_id) > 0
                                 || m_in_opposite_segments.count(turn.operations[1].seg_id) > 0)
                         {
-				add_angles(index, turn.point, turn.operations[0]);
-				add_angles(index, turn.point, turn.operations[1]);
+				add_angles(index);
                         }
                 }
 
+        // Pass 2: add multi intersection
+        add_multi_intersections_to_occupation_map();
+
+        // Pass 3: fill all segments intersecting in points present in the map
         index = 0;
         for (typename boost::range_iterator<turn_vector_type>::type it =
             boost::begin(m_turns); it != boost::end(m_turns); ++it, ++index)
@@ -445,9 +987,7 @@
                 // See if it is in the map
                 if (m_occupation_map.contains(turn.point))
                 {
-//std::cout << "Adding point " << geometry::wkt(turn.point) << std::endl;
-				    add_angles(index, turn.point, turn.operations[0]);
-				    add_angles(index, turn.point, turn.operations[1]);
+				    add_angles(index);
                 }
                         }
                 }
@@ -462,7 +1002,7 @@
         for (typename boost::range_iterator<turn_vector_type>::type it =
             boost::begin(m_turns); it != boost::end(m_turns); ++it)
                 {
-            if (it->count_on_occupied == 0)
+            //if (it->count_on_occupied == 0)
             {
                 typename std::vector<piece>::const_iterator pit;
                 for (pit = boost::begin(m_pieces);
@@ -473,18 +1013,38 @@
                 }
             }
         }
+    }
+
+    template <typename Turn>
+    inline bool classify_turn_inside(Turn const& turn) const
+    {
+        return turn.count_within > 0 
+            || turn.count_on_multi > 0
+            || turn.count_on_helper > 0
+			|| turn.count_on_occupied > 0
+            ;
+    }
 
+    inline void classify_inside()
+    {
                 // Set results:
         for (typename boost::range_iterator<turn_vector_type>::type it =
             boost::begin(m_turns); it != boost::end(m_turns); ++it)
         {
-            if (it->count_within > 0 
-                || it->count_on_helper > 0
-				|| it->count_on_occupied > 0
-				)
+            if (classify_turn_inside(*it))
             {
                 it->location = inside_buffer;
             }
+
+#if defined(BOOST_GEOMETRY_COUNT_DOUBLE_UU)
+            else if (it->count_on_uu > 0)
+            {
+                extern int g_count_double_uu;
+                g_count_double_uu++;
+                std::cout << "UU";
+            }
+#endif
+
         }
     }
 
@@ -533,6 +1093,7 @@
                 get_occupation();
         classify_occupied_locations();
         classify_turns();
+        classify_inside();
 
         if (boost::is_same<typename tag_cast<typename tag<Geometry>::type, areal_tag>::type, areal_tag>())
         {
Modified: trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp
==============================================================================
--- trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp	(original)
+++ trunk/boost/geometry/extensions/algorithms/buffer/buffered_piece_collection_with_mapper.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -35,42 +35,45 @@
     template <typename Mapper>
     inline void map_opposite_locations(Mapper& mapper)
     {
-        for (typename boost::range_iterator<clustered_location_type const>::type it =
-            boost::begin(clustered_turn_locations); 
-            it != boost::end(clustered_turn_locations); ++it)
+        for (typename boost::range_iterator<typename occupation_map_type::map_type>::type it =
+            boost::begin(m_occupation_map.map); 
+            it != boost::end(m_occupation_map.map); ++it)
         {
-			mapper.map(it->first, "fill:rgb(0,128,0);", 3);
+			mapper.map(it->first, it->second.occupied() ? "fill:rgb(255,0,255);" : "fill:rgb(0,192,0);", 7);
 
                     std::ostringstream out;
-			out << it->second.angles.size();
+			out << it->second.angles.size() << std::endl;
             for (std::set<int>::const_iterator sit = it->second.turn_indices.begin(); sit != it->second.turn_indices.end(); ++sit)
             {
                 out << "," << *sit;
             }
-    		mapper.text(it->first, out.str(), "fill:rgb(0,128,0);font-family='Arial';font-size:8px", 6, 8);
+    		mapper.text(it->first, out.str(), "fill:rgb(0,0,0);font-family='Arial';font-size:10px", 6, 8);
 
             for (unsigned int i = 0; i < it->second.angles.size(); i++)
             {
+                double const angle = it->second.angles[i].angle;
+                bool const incoming = it->second.angles[i].incoming;
+                segment_identifier seg_id = it->second.angles[i].seg_id;
+
                 geometry::model::linestring<point_type> line;
-                angle_info const& tp = it->second.angles[i];
                 point_type p1, p2;
-                geometry::set<0>(p1, geometry::get<0>(it->first) + cos(tp.angle) * 0.1);
-                geometry::set<1>(p1, geometry::get<1>(it->first) + sin(tp.angle) * 0.1);
-                geometry::set<0>(p2, geometry::get<0>(it->first) + cos(tp.angle) * 0.4);
-                geometry::set<1>(p2, geometry::get<1>(it->first) + sin(tp.angle) * 0.4);
+                geometry::set<0>(p1, geometry::get<0>(it->first) + cos(angle) * 0.1);
+                geometry::set<1>(p1, geometry::get<1>(it->first) + sin(angle) * 0.1);
+                geometry::set<0>(p2, geometry::get<0>(it->first) + cos(angle) * 0.4);
+                geometry::set<1>(p2, geometry::get<1>(it->first) + sin(angle) * 0.4);
                             std::ostringstream out;
-                //out << tp.angle << " " << int(tp.incoming);
-                out << (tp.incoming ? "i" : "o") << " " << i;
-                if (tp.incoming)
+                out << (incoming ? "i" : "o") << " " << si(seg_id);
+                // out << " " << angle;
+                if (incoming)
                 {
                     int offset = 7;
-                    if (tp.debug) offset += 5;
-                    out << " " << tp.debug_info;
+                    //if (tp.debug) offset += 5;
+                    //out << " " << tp.debug_info;
                     line.push_back(p1);
                     line.push_back(p2);
                                 mapper.map(line, "stroke:rgb(0,0,255);stroke-width:1", 1);
                                 mapper.map(p1, "fill:rgb(0,0,0);", 2);
-    			    mapper.text(p2, out.str(), "fill:rgb(0,0,255);font-family='Arial';font-size:8px", 2, offset);
+    			    mapper.text(p2, out.str(), "fill:rgb(0,0,0);font-family='Arial';font-size:8px", 2, offset);
                 }
                 else
                 {
@@ -78,7 +81,7 @@
                     line.push_back(p2);
                                 mapper.map(line, "stroke:rgb(255,0,0);stroke-width:1", 1);
                                 mapper.map(p2, "fill:rgb(0,0,0);", 2);
-    			    mapper.text(p2, out.str(), "fill:rgb(0,0,255);font-family='Arial';font-size:8px", 2, -2);
+    			    mapper.text(p2, out.str(), "fill:rgb(0,0,0);font-family='Arial';font-size:8px", 2, -2);
                 }
             }
         }
@@ -108,7 +111,9 @@
                                         case inside_original : fill = "fill:rgb(0,0,255);"; color = 'b'; break;
                                 }
                                 std::ostringstream out;
-				out << it->operations[0].piece_index << "/" << it->operations[1].piece_index << std::endl;
+				out << it->operations[0].piece_index << "/" << it->operations[1].piece_index 
+                    << " " << si(it->operations[0].seg_id) << "/" << si(it->operations[1].seg_id)
+                    << std::endl;
                                 //out << " " <<  m_pieces[it->operations[0].piece_index].first_seg_id.segment_index
                                 //     << "+" << m_pieces[it->operations[1].piece_index].first_seg_id.segment_index;
                                 //out << " " <<  m_pieces[it->operations[0].piece_index].index
@@ -123,6 +128,7 @@
                                         << "-" << it->count_on_corner
                                         << "-" << it->count_on_offsetted
                                         << "-" << it->count_on_occupied
+					<< "-" << it->count_on_multi
                                         //<< it->debug_string
                                         ;
                                 out << color << std::endl;
Modified: trunk/boost/geometry/extensions/algorithms/buffer/line_line_intersection.hpp
==============================================================================
--- trunk/boost/geometry/extensions/algorithms/buffer/line_line_intersection.hpp	(original)
+++ trunk/boost/geometry/extensions/algorithms/buffer/line_line_intersection.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -44,6 +44,8 @@
 
         coordinate_type denominator = det(x1 - x2, y1 - y2, x3 - x4, y3 - y4);
 
+        // TODO: use something else then denominator (sides?) to determine this.
+
         // If denominator is zero, segments are parallel.
         // We have context information, so know that it should then
         // be the case that line1.p2 == line2.p1, and that is the
Modified: trunk/boost/geometry/extensions/algorithms/buffer/side_on_convex_range.hpp
==============================================================================
--- trunk/boost/geometry/extensions/algorithms/buffer/side_on_convex_range.hpp	(original)
+++ trunk/boost/geometry/extensions/algorithms/buffer/side_on_convex_range.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -22,9 +22,8 @@
 template <int D>
 struct collinear_point_on_segment_check
 {
-    typedef double coordinate_type; // TODO just use one - they're all the same (for buffer)
-
-    static inline bool apply_sorted(coordinate_type const& subject, coordinate_type const& c1, coordinate_type const& c2)
+	template <typename T>
+    static inline bool apply_sorted(T const& subject, T const& c1, T const& c2)
     {
         return subject >= c1 && subject <= c2;
     }
@@ -32,6 +31,7 @@
     template <typename P0, typename P1, typename P2>
     static inline bool apply(P0 const& subject, P1 const& p1, P2 const& p2)
     {
+	    typedef typename geometry::coordinate_type<P0>::type coordinate_type;
         coordinate_type const cs = geometry::get<D>(subject);
         coordinate_type const c1 = geometry::get<D>(p1);
         coordinate_type const c2 = geometry::get<D>(p2);
Modified: trunk/boost/geometry/extensions/strategies/buffer.hpp
==============================================================================
--- trunk/boost/geometry/extensions/strategies/buffer.hpp	(original)
+++ trunk/boost/geometry/extensions/strategies/buffer.hpp	2012-05-26 16:58:56 EDT (Sat, 26 May 2012)
@@ -227,12 +227,13 @@
         dx2 /= buffer_distance;
         dy2 /= buffer_distance;
 
-        promoted_type angle_diff = std::acos(dx1 * dx2 + dy1 * dy2);
+        promoted_type angle_diff = acos(dx1 * dx2 + dy1 * dy2);
 
         // Default might be 100 steps for a full circle (2 pi)
         promoted_type const steps_per_circle = 100.0;
-        int n = int(steps_per_circle * angle_diff 
-                    / (2.0 * geometry::math::pi<promoted_type>()));
+		promoted_type two = 2.0;
+        int n = boost::numeric_cast<int>(steps_per_circle * angle_diff 
+                    / (two * geometry::math::pi<promoted_type>()));
 
                 if (n > 1000)
                 {
@@ -245,7 +246,7 @@
                         return;
                 }
 
-        promoted_type const angle1 = std::atan2(dy1, dx1);
+        promoted_type const angle1 = atan2(dy1, dx1);
         promoted_type diff = angle_diff / promoted_type(n);
         promoted_type a = angle1 - diff;