[geometry] buffer, extracted occupation info to separate file
[boost:cmake.git] / boost / geometry / extensions / algorithms / buffer / buffered_piece_collection.hpp
1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2
3 // Copyright (c) 2012 Barend Gehrels, Amsterdam, the Netherlands.
4
5 // Use, modification and distribution is subject to the Boost Software License,
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
8
9 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_BUFFERED_PIECE_COLLECTION_HPP
10 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_BUFFERED_PIECE_COLLECTION_HPP
11
12
13 #include <algorithm>
14 #include <cstddef>
15 #include <set>
16 #include <boost/range.hpp>
17
18 #include <boost/geometry/core/coordinate_type.hpp>
19 #include <boost/geometry/core/point_type.hpp>
20
21 #include <boost/geometry/algorithms/equals.hpp>
22 #include <boost/geometry/algorithms/covered_by.hpp>
23 #include <boost/geometry/extensions/strategies/buffer_side.hpp>
24
25 #include <boost/geometry/algorithms/detail/overlay/calculate_distance_policy.hpp>
26 #include <boost/geometry/algorithms/detail/overlay/enrichment_info.hpp>
27 #include <boost/geometry/algorithms/detail/overlay/traversal_info.hpp>
28 #include <boost/geometry/algorithms/detail/overlay/traverse.hpp>
29 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
30 #include <boost/geometry/algorithms/detail/occupation_info.hpp>
31 #include <boost/geometry/algorithms/detail/partition.hpp>
32
33
34 #include <boost/geometry/extensions/algorithms/buffer/buffered_ring.hpp>
35 #include <boost/geometry/extensions/algorithms/buffer/buffer_policies.hpp>
36 #include <boost/geometry/extensions/algorithms/buffer/side_on_convex_range.hpp>
37
38
39 namespace boost { namespace geometry
40 {
41
42
43 #ifndef DOXYGEN_NO_DETAIL
44 namespace detail { namespace buffer
45 {
46
47 enum piece_type
48 {
49     buffered_segment, buffered_join, buffered_flat_end
50 };
51
52 enum segment_relation_code
53 {
54     segment_relation_on_left, 
55     segment_relation_on_right, 
56     segment_relation_within,
57     segment_relation_disjoint
58 };
59
60
61 // In the end this will go (if we have a multi-point within/covered_by geometry)
62 // which is optimized for multi-points and skips linestrings
63 template <typename tag>
64 struct check_original
65 {
66 };
67
68 template <>
69 struct check_original<polygon_tag>
70 {
71     template <typename Point, typename Geometry>
72     static inline int apply(Point const& point, Geometry const& geometry)
73     {
74         return geometry::covered_by(point, geometry) ? 1 : -1;
75     }
76 };
77
78 template <>
79 struct check_original<linestring_tag>
80 {
81     template <typename Point, typename Geometry>
82     static inline int apply(Point const& point, Geometry const& geometry)
83     {
84         return 0;
85     }
86 };
87
88
89 template <typename Ring>
90 struct buffered_piece_collection
91 {
92     typedef typename geometry::point_type<Ring>::type point_type;
93     typedef typename geometry::coordinate_type<Ring>::type coordinate_type;
94
95     struct piece
96     {
97         piece_type type;
98         int index;
99
100         // These both form a complete clockwise ring for each piece (with one dupped point)
101
102                 // 1: half, part of offsetted_rings
103         segment_identifier first_seg_id;
104         int last_segment_index; // no segment-identifier - it is always the same
105
106                 // 2: half, not part (will be indexed in one vector too)
107         std::vector<point_type> helper_segments; // 3 points for segment, 2 points for join - 0 points for flat-end
108     };
109
110
111     typedef typename strategy::side::services::default_strategy
112                 <
113                         typename cs_tag<point_type>::type
114                 >::type side_strategy;
115     typedef std::vector<piece> piece_vector;
116
117     piece_vector m_pieces;
118     buffered_ring_collection<buffered_ring<Ring> > offsetted_rings; // indexed by multi_index
119     buffered_ring_collection<Ring> traversed_rings;
120     segment_identifier current_segment_id;
121
122     typedef std::vector<buffer_turn_info<point_type> > turn_vector_type;
123     typedef detail::overlay::get_turn_info
124         <
125             point_type, point_type, buffer_turn_info<point_type>,
126             turn_assign_for_buffer
127         > turn_policy;
128     turn_vector_type m_turns;
129
130
131         // To check clustered locations we keep track of segments being opposite somewhere
132         std::set<segment_identifier> m_in_opposite_segments;
133
134         struct buffer_occupation_info : public occupation_info<angle_info<coordinate_type> >
135         {
136                 std::set<segment_identifier> seg_ids;
137                 std::set<int> turn_indices;
138         };
139
140     typedef occupation_map<point_type, buffer_occupation_info> occupation_map_type;
141     occupation_map_type m_occupation_map;
142
143
144         struct redundant_turn
145         {
146                 inline bool operator()(buffer_turn_info<point_type> const& turn) const
147                 {
148                         // Erase discarded turns (location not OK) and the turns
149                         // only used to detect oppositeness.
150                         return turn.location != location_ok 
151                                 || turn.opposite();
152                 }
153         };
154
155
156     inline bool is_neighbor(piece const& piece1, piece const& piece2) const
157     {
158         if (std::abs(piece1.index - piece2.index) == 1)
159         {
160             return true;
161         }
162
163         // TODO: of the same multi_index
164         int const b = boost::size(m_pieces) - 1; // back
165         return (piece1.index == 0 && piece2.index == b)
166             || (piece1.index == b && piece2.index == 0)
167             ;
168     }
169
170     inline bool skip_neighbor(piece const& piece1, piece const& piece2) const
171     {
172         return piece1.type != piece2.type && is_neighbor(piece1, piece2);
173     }
174
175     template <typename Range, typename Iterator>
176     inline typename boost::range_value<Range const>::type next_point(Range const& range, Iterator it) const
177     {
178         Iterator next = it;
179         ++next;
180
181                 // Get next point. If end get second point (because ring is assumed to be closed)
182                 return next != boost::end(range) ? *next : *(boost::begin(range) + 1);
183     }
184
185     inline void calculate_turns(piece const& piece1, piece const& piece2)
186     {
187         buffer_turn_info<point_type> the_model;
188         the_model.operations[0].piece_index = piece1.index;
189         the_model.operations[0].seg_id = piece1.first_seg_id;
190
191         typedef typename boost::range_iterator<buffered_ring<Ring> const>::type iterator;
192
193                 segment_identifier seg_id1 = piece1.first_seg_id;
194                 segment_identifier seg_id2 = piece2.first_seg_id;
195
196                 if (seg_id1.segment_index < 0 || seg_id2.segment_index < 0)
197                 {
198                         return;
199                 }
200
201                 buffered_ring<Ring> const& ring1 = offsetted_rings[seg_id1.multi_index];
202                 iterator it1_first = boost::begin(ring1) + seg_id1.segment_index;
203                 iterator it1_last = boost::begin(ring1) + piece1.last_segment_index;
204
205                 buffered_ring<Ring> const& ring2 = offsetted_rings[seg_id2.multi_index];
206                 iterator it2_first = boost::begin(ring2) + seg_id2.segment_index;
207                 iterator it2_last = boost::begin(ring2) + piece2.last_segment_index;
208
209         iterator it1 = it1_first;
210         for (iterator prev1 = it1++; 
211                 it1 != it1_last; 
212                 prev1 = it1++, the_model.operations[0].seg_id.segment_index++)
213         {
214             the_model.operations[1].piece_index = piece2.index;
215             the_model.operations[1].seg_id = piece2.first_seg_id;
216
217             iterator it2 = it2_first;
218             for (iterator prev2 = it2++; 
219                     it2 != it2_last;
220                     prev2 = it2++, the_model.operations[1].seg_id.segment_index++)
221             {
222                 // Revert (this is used more often - should be common function TODO)
223                 the_model.operations[0].other_id = the_model.operations[1].seg_id;
224                 the_model.operations[1].other_id = the_model.operations[0].seg_id;
225
226                 turn_policy::apply(*prev1, *it1, next_point(ring1, it1),
227                                     *prev2, *it2, next_point(ring2, it2),
228                                     the_model, std::back_inserter(m_turns));
229             }
230         }
231     }
232
233         inline void fill_opposite_segments()
234         {
235         for (typename boost::range_iterator<turn_vector_type const>::type it =
236             boost::begin(m_turns); it != boost::end(m_turns); ++it)
237         {
238                         if (it->is_opposite)
239                         {
240                                 m_in_opposite_segments.insert(it->operations[0].seg_id);
241                                 m_in_opposite_segments.insert(it->operations[1].seg_id);
242                         }
243                 }
244         }
245
246         inline segment_relation_code get_segment_relation(point_type const& point,
247                 segment_identifier const& seg_id) const
248         {
249                 typedef typename boost::range_iterator<std::vector<point_type> const>::type iterator_type;
250                 iterator_type it = boost::begin(offsetted_rings[seg_id.multi_index]) + seg_id.segment_index;
251                 iterator_type prev = it++;
252                 int side = side_strategy::apply(point, *prev, *it);
253         if (side == 0)
254         {
255             if (geometry::equals(point, *prev))
256             {
257                 return segment_relation_on_left;
258             }
259             else if (geometry::equals(point, *it))
260             {
261                 return segment_relation_on_right;
262             }
263             else if (collinear_point_on_segment(point, *prev, *it))
264             {
265                 return segment_relation_within;
266             }
267         }
268         return segment_relation_disjoint;
269         }
270
271     inline void add_angles(int index, point_type const& point, buffer_turn_operation<point_type> const& operation)
272     {
273         buffer_occupation_info& info = m_occupation_map.find_or_insert(point);
274         info.turn_indices.insert(index);
275         if (info.seg_ids.count(operation.seg_id) <= 0)
276         {
277             info.seg_ids.insert(operation.seg_id);
278             add_incoming_and_outgoing_angles(point, 
279                                                 offsetted_rings[operation.seg_id.multi_index], 
280                                                 operation.seg_id.segment_index, 
281                                                 info);
282         }
283     }
284
285     inline void classify_turn(buffer_turn_info<point_type>& turn, piece const& pc) const
286     {
287         if (pc.type == buffered_flat_end)
288         {
289             return;
290         }
291
292                 int flat_ends_involved = 0;
293         for (int i = 0; i < boost::size(turn.operations); i++)
294         {
295                         // Don't check any turn against a piece of which is itself the result
296                         if (turn.operations[i].piece_index == pc.index)
297                         {
298                                 return;
299                         }
300
301             piece const& piece_from_intersection = m_pieces[turn.operations[i].piece_index];
302             if (piece_from_intersection.type == buffered_flat_end)
303                         {
304                                 flat_ends_involved++;
305                         }
306                 }
307
308         int side_helper = side_on_convex_range<side_strategy>(turn.point, pc.helper_segments);
309         if (side_helper == 1)
310         {
311                         // Left or outside
312             return;
313         }
314
315                 segment_identifier seg_id = pc.first_seg_id;
316                 if (seg_id.segment_index < 0)
317                 {
318                         // Should not occur
319                         std::cout << "Warning: negative segment_index" << std::endl;
320                         return;
321                 }
322
323                 segment_identifier on_segment_seg_id;
324
325                 buffered_ring<Ring> const& ring = offsetted_rings[seg_id.multi_index];
326
327         int const side_offsetted = side_on_convex_range<side_strategy>(turn.point, 
328                                                 boost::begin(ring) + seg_id.segment_index, 
329                                                 boost::begin(ring) + pc.last_segment_index,
330                                                 seg_id, on_segment_seg_id);
331         if (side_offsetted == 1)
332         {
333             return;
334         }
335
336         if (side_offsetted == -1 && side_helper == -1)
337         {
338             // It is within (assumed that both halves form a closed convex clockwise ring)
339             turn.count_within++;
340         }
341         if (side_offsetted == 0)
342         {
343                 turn.count_on_offsetted++;
344         }
345         if (side_helper == 0)
346         {
347             if (geometry::equals(turn.point, pc.helper_segments.back())
348                 || geometry::equals(turn.point, pc.helper_segments.front()))
349             {
350                 turn.count_on_corner++;
351             }
352             else
353             {
354                                 if (flat_ends_involved == 0)
355                                 {
356                                         turn.count_on_helper++;
357 #ifdef BOOST_GEOMETRY_DEBUG_WITH_MAPPER
358                                         std::ostringstream out;
359                                         out << "HLP " << pc.index;
360                                         turn.debug_string += out.str();
361 #endif
362                                 }
363                                 else
364                                 {
365                                         turn.count_on_corner++;
366                                 }
367             }
368         }
369     }
370
371         inline void classify_occupied_locations()
372         {
373         for (typename boost::range_iterator<typename occupation_map_type::map_type>::type it =
374                     boost::begin(m_occupation_map.map);
375                         it != boost::end(m_occupation_map.map); ++it)
376         {
377                 buffer_occupation_info& info = it->second;
378                         if (info.occupied())
379                         {
380                                 for (std::set<int>::const_iterator sit = info.turn_indices.begin();
381                                         sit != info.turn_indices.end();
382                                         ++sit)
383                                 {
384                                         m_turns[*sit].count_on_occupied++;
385                                 }
386                         }
387         }
388         }
389
390         inline void get_occupation()
391     {
392                 fill_opposite_segments();
393
394         int index = 0;
395         for (typename boost::range_iterator<turn_vector_type>::type it =
396             boost::begin(m_turns); it != boost::end(m_turns); ++it, ++index)
397         {
398                         buffer_turn_info<point_type>& turn = *it;
399                         if (m_in_opposite_segments.count(turn.operations[0].seg_id) > 0
400                                 || m_in_opposite_segments.count(turn.operations[1].seg_id) > 0)
401                         {
402                                 add_angles(index, turn.point, turn.operations[0]);
403                                 add_angles(index, turn.point, turn.operations[1]);
404                         }
405                 }
406         }
407
408     inline void classify_turns()
409     {
410
411         // Check if it is inside any of the pieces
412         // Now: quadratic
413         // TODO: in partition.
414         for (typename boost::range_iterator<turn_vector_type>::type it =
415             boost::begin(m_turns); it != boost::end(m_turns); ++it)
416                 {
417             if (it->count_on_occupied == 0)
418             {
419                 typename std::vector<piece>::const_iterator pit;
420                 for (pit = boost::begin(m_pieces);
421                     pit != boost::end(m_pieces);
422                     ++pit)
423                 {
424                     classify_turn(*it, *pit);
425                 }
426             }
427         }
428
429                 // Set results:
430         for (typename boost::range_iterator<turn_vector_type>::type it =
431             boost::begin(m_turns); it != boost::end(m_turns); ++it)
432         {
433             if (it->count_within > 0 
434                 || it->count_on_helper > 0
435                                 || it->count_on_occupied > 0
436                                 )
437             {
438                 it->location = inside_buffer;
439             }
440         }
441     }
442
443     template <typename Geometry>
444     inline void check_remaining_points(Geometry const& input_geometry, int factor)
445     {
446         // TODO: this should be done as a collection-of-points, for performance
447         for (typename boost::range_iterator<turn_vector_type>::type it =
448             boost::begin(m_turns); it != boost::end(m_turns); ++it)
449         {
450             if (it->location == location_ok)
451                         {
452                                 int code = check_original
453                         <
454                             typename geometry::tag<Geometry>::type
455                         >::apply(it->point, input_geometry);
456                                 if (code * factor == 1)
457                                 {
458                                         it->location = inside_original;
459                                 }
460                         }
461         }
462     }
463
464     template <typename Geometry>
465     inline void get_turns(Geometry const& input_geometry, int factor)
466     {
467                 // Now: quadratic
468         // TODO use partition
469
470         for(typename piece_vector::const_iterator it1 = boost::begin(m_pieces);
471             it1 != boost::end(m_pieces);
472             ++it1)
473         {
474             for(typename piece_vector::const_iterator it2 = it1 + 1;
475                 it2 != boost::end(m_pieces);
476                 ++it2)
477             {
478                 if (! skip_neighbor(*it1, *it2))
479                 {
480                     calculate_turns(*it1, *it2);
481                 }
482             }
483         }
484
485                 get_occupation();
486         classify_occupied_locations();
487         classify_turns();
488
489         if (boost::is_same<typename tag_cast<typename tag<Geometry>::type, areal_tag>::type, areal_tag>())
490         {
491             check_remaining_points(input_geometry, factor);
492         }
493     }
494
495     inline void start_new_ring()
496     {
497         int const n = offsetted_rings.size();
498         current_segment_id.source_index = 0;
499         current_segment_id.multi_index = n;
500         current_segment_id.ring_index = -1;
501         current_segment_id.segment_index = 0;
502
503         offsetted_rings.resize(n + 1);
504     }
505
506     inline int add_point(point_type const& p)
507     {
508         BOOST_ASSERT
509             (
510                 boost::size(offsetted_rings) > 0
511             );
512
513         current_segment_id.segment_index++;
514         offsetted_rings.back().push_back(p);
515                 return offsetted_rings.back().size();
516     }
517
518     inline piece& add_piece(piece_type type, bool decrease_by_one)
519     {
520         piece pc;
521         pc.type = type;
522         pc.index = boost::size(m_pieces);
523         pc.first_seg_id = current_segment_id;
524
525         std::size_t const n = boost::size(offsetted_rings.back());
526         pc.first_seg_id.segment_index = decrease_by_one ? n - 1 : n;
527
528         m_pieces.push_back(pc);
529         return m_pieces.back();
530     }
531
532     inline void add_piece(piece_type type, point_type const& p1, point_type const& p2, 
533             point_type const& b1, point_type const& b2)
534     {
535         bool const last_type_join = ! m_pieces.empty() 
536                                 && m_pieces.back().first_seg_id.multi_index == current_segment_id.multi_index
537                 && m_pieces.back().type == buffered_join;
538
539         piece& pc = add_piece(type, last_type_join);
540
541         // If it follows the same piece-type point both should be added.
542         // There should be two intersections later and it should be discarded.
543         // But for need it to calculate intersections
544         if (! last_type_join)
545         {
546             add_point(b1);
547         }
548         pc.last_segment_index = add_point(b2);
549
550         pc.helper_segments.push_back(b2);
551         pc.helper_segments.push_back(p2);
552         pc.helper_segments.push_back(p1);
553         pc.helper_segments.push_back(b1);
554     }
555
556     template <typename Range>
557     inline piece& add_piece(piece_type type, Range const& range)
558     {
559         piece& pc = add_piece(type, true);
560
561         bool first = true;
562                 int last = offsetted_rings.back().size() + 1;
563         for (typename Range::const_iterator it = boost::begin(range);
564             it != boost::end(range);
565             ++it)
566         {
567             bool add = true;
568             if (first)
569             {
570                 // Only for very first one, add first. In all other cases it is shared with previous.
571                 add = offsetted_rings.back().empty();
572                 first = false;
573             }
574             if (add)
575             {
576                 last = add_point(*it);
577             }
578
579         }
580
581         pc.last_segment_index = last;
582
583         return pc;
584     }
585
586     template <typename Range>
587     inline void add_piece(piece_type type, point_type const& p, Range const& range)
588     {
589         piece& pc = add_piece(type, range);
590
591         if (boost::size(range) > 0)
592         {
593             pc.helper_segments.push_back(range.back());
594             pc.helper_segments.push_back(p);
595             pc.helper_segments.push_back(range.front());
596         }
597     }
598
599     inline void enrich()
600     {
601         typedef typename strategy::side::services::default_strategy
602         <
603             typename cs_tag<Ring>::type
604         >::type side_strategy_type;
605
606         enrich_intersection_points<false, false>(m_turns,
607                     detail::overlay::operation_union,
608                     offsetted_rings, offsetted_rings,
609                     side_strategy_type());
610     }
611
612     // Discards all rings which do have not-OK intersection points only.
613     // Those can never be traversed and should not be part of the output.
614     inline void discard_rings()
615     {
616         for (typename boost::range_iterator<turn_vector_type const>::type it =
617             boost::begin(m_turns); it != boost::end(m_turns); ++it)
618         {
619             if (it->location != location_ok)
620             {
621                 offsetted_rings[it->operations[0].seg_id.multi_index].has_discarded_intersections = true;
622                 offsetted_rings[it->operations[1].seg_id.multi_index].has_discarded_intersections = true;
623             }
624             else
625             {
626                 offsetted_rings[it->operations[0].seg_id.multi_index].has_accepted_intersections = true;
627                 offsetted_rings[it->operations[1].seg_id.multi_index].has_accepted_intersections = true;
628             }
629         }
630     }
631                     
632     inline void discard_turns()
633     {
634         m_turns.erase
635                         (
636                                 std::remove_if(boost::begin(m_turns), boost::end(m_turns),
637                                                                 redundant_turn()),
638                                 boost::end(m_turns)
639                         );
640
641     }
642
643     inline void traverse()
644     {
645         typedef detail::overlay::traverse
646             <
647                 false, false, 
648                 buffered_ring_collection<buffered_ring<Ring> >, buffered_ring_collection<buffered_ring<Ring > >,
649                 backtrack_for_buffer
650             > traverser;
651
652         traversed_rings.clear();
653         traverser::apply(offsetted_rings, offsetted_rings,
654                         detail::overlay::operation_union,
655                         m_turns, traversed_rings);
656     }
657
658     template <typename GeometryOutput, typename OutputIterator>
659     inline OutputIterator assign(OutputIterator out)
660     {
661         typedef detail::overlay::ring_properties<point_type> properties;
662
663         std::map<ring_identifier, properties> selected;
664
665         // Select all rings which do not have any self-intersection (other ones should be traversed)
666         int index = 0;
667         for(typename buffered_ring_collection<buffered_ring<Ring> >::const_iterator it = boost::begin(offsetted_rings);
668             it != boost::end(offsetted_rings);
669             ++it, ++index)
670         {
671             if (! it->has_intersections())
672             {
673                 ring_identifier id(0, index, -1);
674                 selected[id] = properties(*it, true);
675             }
676         }
677
678         // Select all created rings
679         index = 0;
680         for (typename boost::range_iterator<buffered_ring_collection<Ring> const>::type
681                 it = boost::begin(traversed_rings);
682                 it != boost::end(traversed_rings);
683                 ++it, ++index)
684         {
685             ring_identifier id(2, index, -1);
686             selected[id] = properties(*it, true);
687         }
688
689         detail::overlay::assign_parents(offsetted_rings, traversed_rings, selected);
690         return detail::overlay::add_rings<GeometryOutput>(selected, offsetted_rings, traversed_rings, out);
691     }
692
693 };
694
695
696 }} // namespace detail::buffer
697 #endif // DOXYGEN_NO_DETAIL
698
699
700 }} // namespace boost::geometry
701
702 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_BUFFERED_PIECE_COLLECTION_HPP