boost/geometry/algorithms/detail/sections/sectionalize.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4 // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
   5 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
   6 // Copyright (c) 2014-2015 Adam Wulkiewicz, Lodz, Poland.
   7
   8 // This file was modified by Oracle on 2013, 2014, 2015, 2017, 2018, 2019.
   9 // Modifications copyright (c) 2013-2019 Oracle and/or its affiliates.
  10
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  12 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
  13
  14 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  15 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  16
  17 // Use, modification and distribution is subject to the Boost Software License,
  18 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  19 // http://www.boost.org/LICENSE_1_0.txt)
  20
  21 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_SECTIONS_SECTIONALIZE_HPP
  22 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_SECTIONS_SECTIONALIZE_HPP
  23
  24 #include <cstddef>
  25 #include <vector>
  26
  27 #include <boost/concept/requires.hpp>
  28 #include <boost/core/ignore_unused.hpp>
  29 #include <boost/mpl/assert.hpp>
  30 #include <boost/mpl/vector_c.hpp>
  31 #include <boost/range.hpp>
  32 #include <boost/static_assert.hpp>
  33 #include <boost/type_traits/is_same.hpp>
  34 #include <boost/type_traits/is_fundamental.hpp>
  35
  36 #include <boost/geometry/core/config.hpp>
  37
  38 #include <boost/geometry/algorithms/assign.hpp>
  39 #include <boost/geometry/algorithms/envelope.hpp>
  40 #include <boost/geometry/algorithms/expand.hpp>
  41
  42 #include <boost/geometry/algorithms/detail/interior_iterator.hpp>
  43 #include <boost/geometry/algorithms/detail/recalculate.hpp>
  44 #include <boost/geometry/algorithms/detail/ring_identifier.hpp>
  45 #include <boost/geometry/algorithms/detail/signed_size_type.hpp>
  46
  47 #include <boost/geometry/core/access.hpp>
  48 #include <boost/geometry/core/closure.hpp>
  49 #include <boost/geometry/core/exterior_ring.hpp>
  50 #include <boost/geometry/core/point_order.hpp>
  51 #include <boost/geometry/core/tags.hpp>
  52
  53 #include <boost/geometry/geometries/concepts/check.hpp>
  54 #include <boost/geometry/util/math.hpp>
  55 #include <boost/geometry/policies/robustness/no_rescale_policy.hpp>
  56 #include <boost/geometry/policies/robustness/robust_point_type.hpp>
  57 #include <boost/geometry/views/closeable_view.hpp>
  58 #include <boost/geometry/views/reversible_view.hpp>
  59 #include <boost/geometry/geometries/segment.hpp>
  60
  61 #include <boost/geometry/algorithms/detail/expand_by_epsilon.hpp>
  62 #include <boost/geometry/algorithms/detail/buffer/buffer_box.hpp>
  63
  64 #include <boost/geometry/strategies/envelope.hpp>
  65 #include <boost/geometry/strategies/expand.hpp>
  66
  67 namespace boost { namespace geometry
  68 {
  69
  70
  71 /*!
  72     \brief Structure containing section information
  73     \details Section information consists of a bounding box, direction
  74         information (if it is increasing or decreasing, per dimension),
  75         index information (begin-end, ring, multi) and the number of
  76         segments in this section
  77
  78     \tparam Box box-type
  79     \tparam DimensionCount number of dimensions for this section
  80     \ingroup sectionalize
  81  */
  82 template
  83 <
  84     typename Box,
  85     std::size_t DimensionCount
  86 >
  87 struct section
  88 {
  89     typedef Box box_type;
  90     static std::size_t const dimension_count = DimensionCount;
  91
  92     int directions[DimensionCount];
  93     ring_identifier ring_id;
  94     Box bounding_box;
  95
  96     // NOTE: size_type could be passed as template parameter
  97     // NOTE: these probably also could be of type std::size_t
  98     signed_size_type begin_index;
  99     signed_size_type end_index;
 100     std::size_t count;
 101     std::size_t range_count;
 102     bool duplicate;
 103     signed_size_type non_duplicate_index;
 104
 105     bool is_non_duplicate_first;
 106     bool is_non_duplicate_last;
 107
 108     inline section()
 109         : begin_index(-1)
 110         , end_index(-1)
 111         , count(0)
 112         , range_count(0)
 113         , duplicate(false)
 114         , non_duplicate_index(-1)
 115         , is_non_duplicate_first(false)
 116         , is_non_duplicate_last(false)
 117     {
 118         assign_inverse(bounding_box);
 119         for (std::size_t i = 0; i < DimensionCount; i++)
 120         {
 121             directions[i] = 0;
 122         }
 123     }
 124 };
 125
 126
 127 /*!
 128     \brief Structure containing a collection of sections
 129     \note Derived from a vector, proves to be faster than of deque
 130     \note vector might be templated in the future
 131     \ingroup sectionalize
 132  */
 133 template <typename Box, std::size_t DimensionCount>
 134 struct sections : std::vector<section<Box, DimensionCount> >
 135 {
 136     typedef Box box_type;
 137     static std::size_t const value = DimensionCount;
 138 };
 139
 140
 141 #ifndef DOXYGEN_NO_DETAIL
 142 namespace detail { namespace sectionalize
 143 {
 144
 145 // NOTE: This utility will NOT work for latitudes, dimension 1 in spherical
 146 // and geographic coordinate system because in these coordinate systems
 147 // e.g. a segment on northern hemisphere may go towards greater latitude
 148 // and then towards lesser latitude.
 149 template
 150 <
 151     typename Point,
 152     typename DimensionVector,
 153     std::size_t Index,
 154     std::size_t Count,
 155     typename CastedCSTag = typename tag_cast
 156                             <
 157                                 typename cs_tag<Point>::type,
 158                                 spherical_tag
 159                             >::type
 160 >
 161 struct get_direction_loop
 162 {
 163     typedef typename boost::mpl::at_c<DimensionVector, Index>::type dimension;
 164
 165     template <typename Segment>
 166     static inline void apply(Segment const& seg,
 167                 int directions[Count])
 168     {
 169         typedef typename coordinate_type<Segment>::type coordinate_type;
 170
 171         coordinate_type const c0 = geometry::get<0, dimension::value>(seg);
 172         coordinate_type const c1 = geometry::get<1, dimension::value>(seg);
 173
 174         directions[Index] = c1 > c0 ? 1 : c1 < c0 ? -1 : 0;
 175
 176         get_direction_loop
 177         <
 178             Point,
 179             DimensionVector,
 180             Index + 1,
 181             Count,
 182             CastedCSTag
 183         >::apply(seg, directions);
 184     }
 185 };
 186
 187 template
 188 <
 189     typename Point,
 190     typename DimensionVector,
 191     std::size_t Count
 192 >
 193 struct get_direction_loop<Point, DimensionVector, 0, Count, spherical_tag>
 194 {
 195     typedef typename boost::mpl::at_c<DimensionVector, 0>::type dimension;
 196
 197     template <typename Segment>
 198     static inline void apply(Segment const& seg,
 199                 int directions[Count])
 200     {
 201         typedef typename coordinate_type<Segment>::type coordinate_type;
 202         typedef typename coordinate_system<Point>::type::units units_t;
 203
 204         coordinate_type const diff = math::longitude_distance_signed
 205                                         <
 206                                             units_t, coordinate_type
 207                                         >(geometry::get<0, 0>(seg),
 208                                           geometry::get<1, 0>(seg));
 209
 210         coordinate_type zero = coordinate_type();
 211         directions[0] = diff > zero ? 1 : diff < zero ? -1 : 0;
 212
 213         get_direction_loop
 214         <
 215             Point,
 216             DimensionVector,
 217             1,
 218             Count,
 219             spherical_tag
 220         >::apply(seg, directions);
 221     }
 222 };
 223
 224 template
 225 <
 226     typename Point,
 227     typename DimensionVector,
 228     std::size_t Count,
 229     typename CastedCSTag
 230 >
 231 struct get_direction_loop<Point, DimensionVector, Count, Count, CastedCSTag>
 232 {
 233     template <typename Segment>
 234     static inline void apply(Segment const&, int [Count])
 235     {}
 236 };
 237
 238
 239 //! Copy one static array to another
 240 template <typename T, std::size_t Index, std::size_t Count>
 241 struct copy_loop
 242 {
 243     static inline void apply(T const source[Count], T target[Count])
 244     {
 245         target[Index] = source[Index];
 246         copy_loop<T, Index + 1, Count>::apply(source, target);
 247     }
 248 };
 249
 250 template <typename T, std::size_t Count>
 251 struct copy_loop<T, Count, Count>
 252 {
 253     static inline void apply(T const [Count], T [Count])
 254     {}
 255 };
 256
 257 //! Compare two static arrays
 258 template <typename T, std::size_t Index, std::size_t Count>
 259 struct compare_loop
 260 {
 261     static inline bool apply(T const array1[Count], T const array2[Count])
 262     {
 263         return array1[Index] != array2[Index]
 264             ? false
 265             : compare_loop
 266                 <
 267                     T, Index + 1, Count
 268                 >::apply(array1, array2);
 269     }
 270 };
 271
 272 template <typename T, std::size_t Count>
 273 struct compare_loop<T, Count, Count>
 274 {
 275     static inline bool apply(T const [Count], T const [Count])
 276     {
 277
 278         return true;
 279     }
 280 };
 281
 282
 283 template <std::size_t Dimension, std::size_t DimensionCount>
 284 struct check_duplicate_loop
 285 {
 286     template <typename Segment>
 287     static inline bool apply(Segment const& seg)
 288     {
 289         if (! geometry::math::equals
 290                 (
 291                     geometry::get<0, Dimension>(seg),
 292                     geometry::get<1, Dimension>(seg)
 293                 )
 294             )
 295         {
 296             return false;
 297         }
 298
 299         return check_duplicate_loop
 300         <
 301                 Dimension + 1, DimensionCount
 302         >::apply(seg);
 303     }
 304 };
 305
 306 template <std::size_t DimensionCount>
 307 struct check_duplicate_loop<DimensionCount, DimensionCount>
 308 {
 309     template <typename Segment>
 310     static inline bool apply(Segment const&)
 311     {
 312         return true;
 313     }
 314 };
 315
 316 //! Assign a value to a static array
 317 template <typename T, std::size_t Index, std::size_t Count>
 318 struct assign_loop
 319 {
 320     static inline void apply(T dims[Count], T const value)
 321     {
 322         dims[Index] = value;
 323         assign_loop<T, Index + 1, Count>::apply(dims, value);
 324     }
 325 };
 326
 327 template <typename T, std::size_t Count>
 328 struct assign_loop<T, Count, Count>
 329 {
 330     static inline void apply(T [Count], T const)
 331     {
 332     }
 333 };
 334
 335 template <typename CSTag>
 336 struct box_first_in_section
 337 {
 338     template <typename Box, typename Point, typename EnvelopeStrategy>
 339     static inline void apply(Box & box, Point const& prev, Point const& curr,
 340                              EnvelopeStrategy const& strategy)
 341     {
 342         geometry::model::referring_segment<Point const> seg(prev, curr);
 343         geometry::envelope(seg, box, strategy);
 344     }
 345 };
 346
 347 template <>
 348 struct box_first_in_section<cartesian_tag>
 349 {
 350     template <typename Box, typename Point, typename ExpandStrategy>
 351     static inline void apply(Box & box, Point const& prev, Point const& curr,
 352                              ExpandStrategy const& )
 353     {
 354         geometry::envelope(prev, box);
 355         geometry::expand(box, curr);
 356     }
 357 };
 358
 359 template <typename CSTag>
 360 struct box_next_in_section
 361 {
 362     template <typename Box, typename Point, typename Strategy>
 363     static inline void apply(Box & box, Point const& prev, Point const& curr,
 364                              Strategy const& strategy)
 365     {
 366         geometry::model::referring_segment<Point const> seg(prev, curr);
 367         geometry::expand(box, seg, strategy);
 368     }
 369 };
 370
 371 template <>
 372 struct box_next_in_section<cartesian_tag>
 373 {
 374     template <typename Box, typename Point, typename Strategy>
 375     static inline void apply(Box & box, Point const& , Point const& curr,
 376                              Strategy const& )
 377     {
 378         geometry::expand(box, curr);
 379     }
 380 };
 381
 382 /// @brief Helper class to create sections of a part of a range, on the fly
 383 template
 384 <
 385     typename Point,
 386     typename DimensionVector
 387 >
 388 struct sectionalize_part
 389 {
 390     static const std::size_t dimension_count
 391         = boost::mpl::size<DimensionVector>::value;
 392
 393     template
 394     <
 395         typename Iterator,
 396         typename RobustPolicy,
 397         typename Sections
 398     >
 399     static inline void apply(Sections& sections,
 400                              Iterator begin, Iterator end,
 401                              RobustPolicy const& robust_policy,
 402                              ring_identifier ring_id,
 403                              std::size_t max_count)
 404     {
 405         typedef typename strategy::envelope::services::default_strategy
 406             <
 407                 segment_tag,
 408                 typename cs_tag<typename Sections::box_type>::type
 409             >::type envelope_strategy_type;
 410
 411         typedef typename strategy::expand::services::default_strategy
 412             <
 413                 segment_tag,
 414                 typename cs_tag<typename Sections::box_type>::type
 415             >::type expand_strategy_type;
 416
 417         apply(sections, begin, end,
 418               robust_policy,
 419               envelope_strategy_type(),
 420               expand_strategy_type(),
 421               ring_id, max_count);
 422     }
 423
 424     template
 425     <
 426         typename Iterator,
 427         typename RobustPolicy,
 428         typename Sections,
 429         typename EnvelopeStrategy,
 430         typename ExpandStrategy
 431     >
 432     static inline void apply(Sections& sections,
 433                              Iterator begin, Iterator end,
 434                              RobustPolicy const& robust_policy,
 435                              EnvelopeStrategy const& envelope_strategy,
 436                              ExpandStrategy const& expand_strategy,
 437                              ring_identifier ring_id,
 438                              std::size_t max_count)
 439     {
 440         boost::ignore_unused(robust_policy);
 441
 442         typedef typename boost::range_value<Sections>::type section_type;
 443         BOOST_STATIC_ASSERT
 444             (
 445                 (static_cast<std::size_t>(section_type::dimension_count)
 446                  == static_cast<std::size_t>(boost::mpl::size<DimensionVector>::value))
 447             );
 448
 449         typedef typename geometry::robust_point_type
 450         <
 451             Point,
 452             RobustPolicy
 453         >::type robust_point_type;
 454
 455         std::size_t const count = std::distance(begin, end);
 456         if (count == 0)
 457         {
 458             return;
 459         }
 460
 461         signed_size_type index = 0;
 462         signed_size_type ndi = 0; // non duplicate index
 463         section_type section;
 464
 465         bool mark_first_non_duplicated = true;
 466         std::size_t last_non_duplicate_index = sections.size();
 467
 468         Iterator it = begin;
 469         robust_point_type previous_robust_point;
 470         geometry::recalculate(previous_robust_point, *it, robust_policy);
 471
 472         for(Iterator previous = it++;
 473             it != end;
 474             ++previous, ++it, index++)
 475         {
 476             robust_point_type current_robust_point;
 477             geometry::recalculate(current_robust_point, *it, robust_policy);
 478             model::referring_segment<robust_point_type> robust_segment(
 479                     previous_robust_point, current_robust_point);
 480
 481             int direction_classes[dimension_count] = {0};
 482             get_direction_loop
 483             <
 484                 Point, DimensionVector, 0, dimension_count
 485             >::apply(robust_segment, direction_classes);
 486
 487             // if "dir" == 0 for all point-dimensions, it is duplicate.
 488             // Those sections might be omitted, if wished, lateron
 489             bool duplicate = false;
 490
 491             if (direction_classes[0] == 0)
 492             {
 493                 // Recheck because ALL dimensions should be checked,
 494                 // not only first one.
 495                 // (dimension_count might be < dimension<P>::value)
 496                 if (check_duplicate_loop
 497                     <
 498                         0, geometry::dimension<Point>::type::value
 499                     >::apply(robust_segment)
 500                     )
 501                 {
 502                     duplicate = true;
 503
 504                     // Change direction-info to force new section
 505                     // Note that wo consecutive duplicate segments will generate
 506                     // only one duplicate-section.
 507                     // Actual value is not important as long as it is not -1,0,1
 508                     assign_loop
 509                     <
 510                         int, 0, dimension_count
 511                     >::apply(direction_classes, -99);
 512                 }
 513             }
 514
 515             if (section.count > 0
 516                 && (! compare_loop
 517                         <
 518                             int, 0, dimension_count
 519                         >::apply(direction_classes, section.directions)
 520                     || section.count > max_count)
 521                 )
 522             {
 523                 if (! section.duplicate)
 524                 {
 525                     last_non_duplicate_index = sections.size();
 526                 }
 527
 528                 sections.push_back(section);
 529                 section = section_type();
 530             }
 531
 532             if (section.count == 0)
 533             {
 534                 section.begin_index = index;
 535                 section.ring_id = ring_id;
 536                 section.duplicate = duplicate;
 537                 section.non_duplicate_index = ndi;
 538                 section.range_count = count;
 539
 540                 if (mark_first_non_duplicated && ! duplicate)
 541                 {
 542                     section.is_non_duplicate_first = true;
 543                     mark_first_non_duplicated = false;
 544                 }
 545
 546                 copy_loop
 547                     <
 548                         int, 0, dimension_count
 549                     >::apply(direction_classes, section.directions);
 550
 551                 // In cartesian this is envelope of previous point expanded with current point
 552                 // in non-cartesian this is envelope of a segment
 553                 box_first_in_section<typename cs_tag<robust_point_type>::type>
 554                     ::apply(section.bounding_box, previous_robust_point, current_robust_point, envelope_strategy);
 555             }
 556             else
 557             {
 558                 // In cartesian this is expand with current point
 559                 // in non-cartesian this is expand with a segment
 560                 box_next_in_section<typename cs_tag<robust_point_type>::type>
 561                     ::apply(section.bounding_box, previous_robust_point, current_robust_point, expand_strategy);
 562             }
 563
 564             section.end_index = index + 1;
 565             section.count++;
 566             if (! duplicate)
 567             {
 568                 ndi++;
 569             }
 570             previous_robust_point = current_robust_point;
 571         }
 572
 573         // Add last section if applicable
 574         if (section.count > 0)
 575         {
 576             if (! section.duplicate)
 577             {
 578                 last_non_duplicate_index = sections.size();
 579             }
 580
 581             sections.push_back(section);
 582         }
 583
 584         if (last_non_duplicate_index < sections.size()
 585             && ! sections[last_non_duplicate_index].duplicate)
 586         {
 587             sections[last_non_duplicate_index].is_non_duplicate_last = true;
 588         }
 589     }
 590 };
 591
 592
 593 template
 594 <
 595     closure_selector Closure,
 596     bool Reverse,
 597     typename Point,
 598     typename DimensionVector
 599 >
 600 struct sectionalize_range
 601 {
 602     template
 603     <
 604         typename Range,
 605         typename RobustPolicy,
 606         typename Sections,
 607         typename EnvelopeStrategy,
 608         typename ExpandStrategy
 609     >
 610     static inline void apply(Range const& range,
 611                              RobustPolicy const& robust_policy,
 612                              Sections& sections,
 613                              EnvelopeStrategy const& envelope_strategy,
 614                              ExpandStrategy const& expand_strategy,
 615                              ring_identifier ring_id,
 616                              std::size_t max_count)
 617     {
 618         typedef typename closeable_view<Range const, Closure>::type cview_type;
 619         typedef typename reversible_view
 620         <
 621             cview_type const,
 622             Reverse ? iterate_reverse : iterate_forward
 623         >::type view_type;
 624
 625         cview_type cview(range);
 626         view_type view(cview);
 627
 628         std::size_t const n = boost::size(view);
 629         if (n == 0)
 630         {
 631             // Zero points, no section
 632             return;
 633         }
 634
 635         if (n == 1)
 636         {
 637             // Line with one point ==> no sections
 638             return;
 639         }
 640
 641         sectionalize_part<Point, DimensionVector>::apply(sections,
 642             boost::begin(view), boost::end(view),
 643             robust_policy, envelope_strategy, expand_strategy,
 644             ring_id, max_count);
 645     }
 646 };
 647
 648 template
 649 <
 650     bool Reverse,
 651     typename DimensionVector
 652 >
 653 struct sectionalize_polygon
 654 {
 655     template
 656     <
 657         typename Polygon,
 658         typename RobustPolicy,
 659         typename Sections,
 660         typename EnvelopeStrategy,
 661         typename ExpandStrategy
 662     >
 663     static inline void apply(Polygon const& poly,
 664                 RobustPolicy const& robust_policy,
 665                 Sections& sections,
 666                 EnvelopeStrategy const& envelope_strategy,
 667                 ExpandStrategy const& expand_strategy,
 668                 ring_identifier ring_id,
 669                 std::size_t max_count)
 670     {
 671         typedef typename point_type<Polygon>::type point_type;
 672         typedef sectionalize_range
 673         <
 674                 closure<Polygon>::value, Reverse,
 675                 point_type, DimensionVector
 676         > per_range;
 677
 678         ring_id.ring_index = -1;
 679         per_range::apply(exterior_ring(poly), robust_policy, sections,
 680                          envelope_strategy, expand_strategy, ring_id, max_count);
 681
 682         ring_id.ring_index++;
 683         typename interior_return_type<Polygon const>::type
 684             rings = interior_rings(poly);
 685         for (typename detail::interior_iterator<Polygon const>::type
 686                 it = boost::begin(rings); it != boost::end(rings); ++it, ++ring_id.ring_index)
 687         {
 688             per_range::apply(*it, robust_policy, sections,
 689                              envelope_strategy, expand_strategy, ring_id, max_count);
 690         }
 691     }
 692 };
 693
 694 template <typename DimensionVector>
 695 struct sectionalize_box
 696 {
 697     template
 698     <
 699         typename Box,
 700         typename RobustPolicy,
 701         typename Sections,
 702         typename EnvelopeStrategy,
 703         typename ExpandStrategy
 704     >
 705     static inline void apply(Box const& box,
 706                 RobustPolicy const& robust_policy,
 707                 Sections& sections,
 708                 EnvelopeStrategy const& envelope_strategy,
 709                 ExpandStrategy const& expand_strategy,
 710                 ring_identifier const& ring_id, std::size_t max_count)
 711     {
 712         typedef typename point_type<Box>::type point_type;
 713
 714         assert_dimension<Box, 2>();
 715
 716         // Add all four sides of the 2D-box as separate section.
 717         // Easiest is to convert it to a polygon.
 718         // However, we don't have the polygon type
 719         // (or polygon would be a helper-type).
 720         // Therefore we mimic a linestring/std::vector of 5 points
 721
 722         // TODO: might be replaced by assign_box_corners_oriented
 723         // or just "convert"
 724         point_type ll, lr, ul, ur;
 725         geometry::detail::assign_box_corners(box, ll, lr, ul, ur);
 726
 727         std::vector<point_type> points;
 728         points.push_back(ll);
 729         points.push_back(ul);
 730         points.push_back(ur);
 731         points.push_back(lr);
 732         points.push_back(ll);
 733
 734         // NOTE: Use cartesian envelope strategy in all coordinate systems
 735         //       because edges of a box are not geodesic segments
 736         sectionalize_range
 737         <
 738                 closed, false,
 739             point_type,
 740             DimensionVector
 741         >::apply(points, robust_policy, sections,
 742                  envelope_strategy, expand_strategy,
 743                  ring_id, max_count);
 744     }
 745 };
 746
 747 template <typename DimensionVector, typename Policy>
 748 struct sectionalize_multi
 749 {
 750     template
 751     <
 752         typename MultiGeometry,
 753         typename RobustPolicy,
 754         typename Sections,
 755         typename EnvelopeStrategy,
 756         typename ExpandStrategy
 757     >
 758     static inline void apply(MultiGeometry const& multi,
 759                 RobustPolicy const& robust_policy,
 760                 Sections& sections,
 761                 EnvelopeStrategy const& envelope_strategy,
 762                 ExpandStrategy const& expand_strategy,
 763                 ring_identifier ring_id,
 764                 std::size_t max_count)
 765     {
 766         ring_id.multi_index = 0;
 767         for (typename boost::range_iterator<MultiGeometry const>::type
 768                     it = boost::begin(multi);
 769             it != boost::end(multi);
 770             ++it, ++ring_id.multi_index)
 771         {
 772             Policy::apply(*it, robust_policy, sections,
 773                           envelope_strategy, expand_strategy,
 774                           ring_id, max_count);
 775         }
 776     }
 777 };
 778
 779 // TODO: If it depends on CS it should probably be made into strategy.
 780 // For now implemented here because of ongoing work on robustness
 781 //   the fact that it interferes with detail::buffer::buffer_box
 782 //   and that we probably need a general strategy for defining epsilon in
 783 //   various coordinate systems, e.g. for point comparison, enlargement of
 784 //   bounding boxes, etc.
 785 template <typename CSTag>
 786 struct expand_by_epsilon
 787     : not_implemented<CSTag>
 788 {};
 789
 790 template <>
 791 struct expand_by_epsilon<cartesian_tag>
 792 {
 793     template <typename Box>
 794     static inline void apply(Box & box)
 795     {
 796         detail::expand_by_epsilon(box);
 797     }
 798 };
 799
 800 template <>
 801 struct expand_by_epsilon<spherical_tag>
 802 {
 803     template <typename Box>
 804     static inline void apply(Box & box)
 805     {
 806         typedef typename coordinate_type<Box>::type coord_t;
 807         static const coord_t eps = boost::is_same<coord_t, float>::value
 808             ? coord_t(1e-6)
 809             : coord_t(1e-12);
 810         detail::expand_by_epsilon(box, eps);
 811     }
 812 };
 813
 814 // TODO: In geographic CS it should probably also depend on FormulaPolicy.
 815 template <>
 816 struct expand_by_epsilon<geographic_tag>
 817     : expand_by_epsilon<spherical_tag>
 818 {};
 819
 820 template <typename Sections, typename Strategy>
 821 inline void enlarge_sections(Sections& sections, Strategy const&)
 822 {
 823     // Enlarge sections slightly, this should be consistent with math::equals()
 824     // and with the tolerances used in general_form intersections.
 825     // This avoids missing turns.
 826
 827     // Points and Segments are equal-compared WRT machine epsilon, but Boxes aren't
 828     // Enlarging Boxes ensures that they correspond to the bound objects,
 829     // Segments in this case, since Sections are collections of Segments.
 830
 831     // It makes section a tiny bit too large, which might cause (a small number)
 832     // of more comparisons
 833     for (typename boost::range_iterator<Sections>::type it = boost::begin(sections);
 834         it != boost::end(sections);
 835         ++it)
 836     {
 837 #if defined(BOOST_GEOMETRY_USE_RESCALING)
 838         detail::sectionalize::expand_by_epsilon
 839             <
 840                 typename Strategy::cs_tag
 841             >::apply(it->bounding_box);
 842
 843 #else
 844         // Expand the box to avoid missing any intersection. The amount is
 845         // should be larger than epsilon. About the value itself: the smaller
 846         // it is, the higher the risk to miss intersections. The larger it is,
 847         // the more comparisons are made. So it should be on the high side.
 848         detail::buffer::buffer_box(it->bounding_box, 0.001, it->bounding_box);
 849 #endif
 850     }
 851 }
 852
 853
 854 }} // namespace detail::sectionalize
 855 #endif // DOXYGEN_NO_DETAIL
 856
 857
 858 #ifndef DOXYGEN_NO_DISPATCH
 859 namespace dispatch
 860 {
 861
 862 template
 863 <
 864     typename Tag,
 865     typename Geometry,
 866     bool Reverse,
 867     typename DimensionVector
 868 >
 869 struct sectionalize
 870 {
 871     BOOST_MPL_ASSERT_MSG
 872         (
 873             false, NOT_OR_NOT_YET_IMPLEMENTED_FOR_THIS_GEOMETRY_TYPE
 874             , (types<Geometry>)
 875         );
 876 };
 877
 878 template
 879 <
 880     typename Box,
 881     bool Reverse,
 882     typename DimensionVector
 883 >
 884 struct sectionalize<box_tag, Box, Reverse, DimensionVector>
 885     : detail::sectionalize::sectionalize_box<DimensionVector>
 886 {};
 887
 888 template
 889 <
 890     typename LineString,
 891     typename DimensionVector
 892 >
 893 struct sectionalize
 894     <
 895         linestring_tag,
 896         LineString,
 897         false,
 898         DimensionVector
 899     >
 900     : detail::sectionalize::sectionalize_range
 901         <
 902             closed, false,
 903             typename point_type<LineString>::type,
 904             DimensionVector
 905         >
 906 {};
 907
 908 template
 909 <
 910     typename Ring,
 911     bool Reverse,
 912     typename DimensionVector
 913 >
 914 struct sectionalize<ring_tag, Ring, Reverse, DimensionVector>
 915     : detail::sectionalize::sectionalize_range
 916         <
 917             geometry::closure<Ring>::value, Reverse,
 918             typename point_type<Ring>::type,
 919             DimensionVector
 920         >
 921 {};
 922
 923 template
 924 <
 925     typename Polygon,
 926     bool Reverse,
 927     typename DimensionVector
 928 >
 929 struct sectionalize<polygon_tag, Polygon, Reverse, DimensionVector>
 930     : detail::sectionalize::sectionalize_polygon
 931         <
 932             Reverse, DimensionVector
 933         >
 934 {};
 935
 936 template
 937 <
 938     typename MultiPolygon,
 939     bool Reverse,
 940     typename DimensionVector
 941 >
 942 struct sectionalize<multi_polygon_tag, MultiPolygon, Reverse, DimensionVector>
 943     : detail::sectionalize::sectionalize_multi
 944         <
 945             DimensionVector,
 946             detail::sectionalize::sectionalize_polygon
 947                 <
 948                     Reverse,
 949                     DimensionVector
 950                 >
 951         >
 952
 953 {};
 954
 955 template
 956 <
 957     typename MultiLinestring,
 958     bool Reverse,
 959     typename DimensionVector
 960 >
 961 struct sectionalize<multi_linestring_tag, MultiLinestring, Reverse, DimensionVector>
 962     : detail::sectionalize::sectionalize_multi
 963         <
 964             DimensionVector,
 965             detail::sectionalize::sectionalize_range
 966                 <
 967                     closed, false,
 968                     typename point_type<MultiLinestring>::type,
 969                     DimensionVector
 970                 >
 971         >
 972
 973 {};
 974
 975 } // namespace dispatch
 976 #endif
 977
 978
 979 /*!
 980     \brief Split a geometry into monotonic sections
 981     \ingroup sectionalize
 982     \tparam Geometry type of geometry to check
 983     \tparam Sections type of sections to create
 984     \param geometry geometry to create sections from
 985     \param robust_policy policy to handle robustness issues
 986     \param sections structure with sections
 987     \param envelope_strategy strategy for envelope calculation
 988     \param expand_strategy strategy for partitions
 989     \param source_index index to assign to the ring_identifiers
 990     \param max_count maximal number of points per section
 991         (defaults to 10, this seems to give the fastest results)
 992
 993  */
 994 template
 995 <
 996     bool Reverse,
 997     typename DimensionVector,
 998     typename Geometry,
 999     typename Sections,
1000     typename RobustPolicy,
1001     typename EnvelopeStrategy,
1002     typename ExpandStrategy
1003 >
1004 inline void sectionalize(Geometry const& geometry,
1005                 RobustPolicy const& robust_policy,
1006                 Sections& sections,
1007                 EnvelopeStrategy const& envelope_strategy,
1008                 ExpandStrategy const& expand_strategy,
1009                 int source_index = 0,
1010                 std::size_t max_count = 10)
1011 {
1012     BOOST_STATIC_ASSERT((! boost::is_fundamental<EnvelopeStrategy>::value));
1013
1014     concepts::check<Geometry const>();
1015
1016     typedef typename boost::range_value<Sections>::type section_type;
1017
1018     // Compiletime check for point type of section boxes
1019     // and point type related to robust policy
1020     typedef typename geometry::coordinate_type
1021     <
1022         typename section_type::box_type
1023     >::type ctype1;
1024     typedef typename geometry::coordinate_type
1025     <
1026         typename geometry::robust_point_type
1027         <
1028             typename geometry::point_type<Geometry>::type,
1029             RobustPolicy
1030         >::type
1031     >::type ctype2;
1032
1033     BOOST_MPL_ASSERT((boost::is_same<ctype1, ctype2>));
1034
1035
1036     sections.clear();
1037
1038     ring_identifier ring_id;
1039     ring_id.source_index = source_index;
1040
1041     dispatch::sectionalize
1042         <
1043             typename tag<Geometry>::type,
1044             Geometry,
1045             Reverse,
1046             DimensionVector
1047         >::apply(geometry, robust_policy, sections,
1048                  envelope_strategy, expand_strategy,
1049                  ring_id, max_count);
1050
1051     detail::sectionalize::enlarge_sections(sections, envelope_strategy);
1052 }
1053
1054
1055 template
1056 <
1057     bool Reverse,
1058     typename DimensionVector,
1059     typename Geometry,
1060     typename Sections,
1061     typename RobustPolicy
1062 >
1063 inline void sectionalize(Geometry const& geometry,
1064                          RobustPolicy const& robust_policy,
1065                          Sections& sections,
1066                          int source_index = 0,
1067                          std::size_t max_count = 10)
1068 {
1069     typedef typename strategy::envelope::services::default_strategy
1070         <
1071             typename tag<Geometry>::type,
1072             typename cs_tag<Geometry>::type
1073         >::type envelope_strategy_type;
1074
1075     typedef typename strategy::expand::services::default_strategy
1076         <
1077             typename boost::mpl::if_c
1078                 <
1079                     boost::is_same<typename tag<Geometry>::type, box_tag>::value,
1080                     box_tag,
1081                     segment_tag
1082                 >::type,
1083             typename cs_tag<Geometry>::type
1084         >::type expand_strategy_type;
1085
1086     boost::geometry::sectionalize
1087         <
1088             Reverse, DimensionVector
1089         >(geometry, robust_policy, sections,
1090           envelope_strategy_type(),
1091           expand_strategy_type(),
1092           source_index, max_count);
1093 }
1094
1095 }} // namespace boost::geometry
1096
1097
1098 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_SECTIONS_SECTIONALIZE_HPP