boost/geometry/algorithms/detail/equals/collect_vectors.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
   4 // Copyright (c) 2008-2014 Bruno Lalande, Paris, France.
   5 // Copyright (c) 2009-2014 Mateusz Loskot, London, UK.
   6 // Copyright (c) 2014-2017 Adam Wulkiewicz, Lodz, Poland.
   7
   8 // This file was modified by Oracle on 2017.
   9 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
  10
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  12
  13 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  14 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  15
  16 // Use, modification and distribution is subject to the Boost Software License,
  17 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  18 // http://www.boost.org/LICENSE_1_0.txt)
  19
  20 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_EQUALS_COLLECT_VECTORS_HPP
  21 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EQUALS_COLLECT_VECTORS_HPP
  22
  23
  24 #include <boost/numeric/conversion/cast.hpp>
  25
  26 #include <boost/geometry/algorithms/detail/interior_iterator.hpp>
  27 #include <boost/geometry/algorithms/detail/normalize.hpp>
  28 #include <boost/geometry/algorithms/not_implemented.hpp>
  29
  30 #include <boost/geometry/core/cs.hpp>
  31 #include <boost/geometry/core/interior_rings.hpp>
  32 #include <boost/geometry/core/tags.hpp>
  33
  34 #include <boost/geometry/formulas/spherical.hpp>
  35
  36 #include <boost/geometry/geometries/concepts/check.hpp>
  37
  38 #include <boost/geometry/util/math.hpp>
  39 #include <boost/geometry/util/range.hpp>
  40
  41 #include <boost/geometry/views/detail/normalized_view.hpp>
  42
  43 #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
  44 #include <boost/geometry/strategies/spherical/ssf.hpp>
  45
  46
  47 namespace boost { namespace geometry
  48 {
  49
  50 // TODO: dispatch only by SideStrategy instead of Geometry/CSTag?
  51
  52 // Since these vectors (though ray would be a better name) are used in the
  53 // implementation of equals() for Areal geometries the internal representation
  54 // should be consistent with the side strategy.
  55 template
  56 <
  57     typename T,
  58     typename Geometry,
  59     typename SideStrategy,
  60     typename CSTag = typename cs_tag<Geometry>::type
  61 >
  62 struct collected_vector
  63     : nyi::not_implemented_tag
  64 {};
  65
  66 // compatible with side_by_triangle cartesian strategy
  67 template <typename T, typename Geometry, typename CT, typename CSTag>
  68 struct collected_vector
  69     <
  70         T, Geometry, strategy::side::side_by_triangle<CT>, CSTag
  71     >
  72 {
  73     typedef T type;
  74
  75     inline collected_vector()
  76     {}
  77
  78     inline collected_vector(T const& px, T const& py,
  79                             T const& pdx, T const& pdy)
  80         : x(px)
  81         , y(py)
  82         , dx(pdx)
  83         , dy(pdy)
  84         //, dx_0(dx)
  85         //, dy_0(dy)
  86     {}
  87
  88     template <typename Point>
  89     inline collected_vector(Point const& p1, Point const& p2)
  90         : x(get<0>(p1))
  91         , y(get<1>(p1))
  92         , dx(get<0>(p2) - x)
  93         , dy(get<1>(p2) - y)
  94         //, dx_0(dx)
  95         //, dy_0(dy)
  96     {}
  97
  98     bool normalize()
  99     {
 100         T magnitude = math::sqrt(
 101             boost::numeric_cast<T>(dx * dx + dy * dy));
 102
 103         // NOTE: shouldn't here math::equals() be called?
 104         if (magnitude > 0)
 105         {
 106             dx /= magnitude;
 107             dy /= magnitude;
 108             return true;
 109         }
 110
 111         return false;
 112     }
 113
 114     // For sorting
 115     inline bool operator<(collected_vector const& other) const
 116     {
 117         if (math::equals(x, other.x))
 118         {
 119             if (math::equals(y, other.y))
 120             {
 121                 if (math::equals(dx, other.dx))
 122                 {
 123                     return dy < other.dy;
 124                 }
 125                 return dx < other.dx;
 126             }
 127             return y < other.y;
 128         }
 129         return x < other.x;
 130     }
 131
 132     inline bool next_is_collinear(collected_vector const& other) const
 133     {
 134         return same_direction(other);
 135     }
 136
 137     // For std::equals
 138     inline bool operator==(collected_vector const& other) const
 139     {
 140         return math::equals(x, other.x)
 141             && math::equals(y, other.y)
 142             && same_direction(other);
 143     }
 144
 145 private:
 146     inline bool same_direction(collected_vector const& other) const
 147     {
 148         // For high precision arithmetic, we have to be
 149         // more relaxed then using ==
 150         // Because 2/sqrt( (0,0)<->(2,2) ) == 1/sqrt( (0,0)<->(1,1) )
 151         // is not always true (at least, it is not for ttmath)
 152         return math::equals_with_epsilon(dx, other.dx)
 153             && math::equals_with_epsilon(dy, other.dy);
 154     }
 155
 156     T x, y;
 157     T dx, dy;
 158     //T dx_0, dy_0;
 159 };
 160
 161 // Compatible with spherical_side_formula which currently
 162 // is the default spherical and geographical strategy
 163 template <typename T, typename Geometry, typename CT, typename CSTag>
 164 struct collected_vector
 165     <
 166         T, Geometry, strategy::side::spherical_side_formula<CT>, CSTag
 167     >
 168 {
 169     typedef T type;
 170
 171     typedef typename coordinate_system<Geometry>::type cs_type;
 172     typedef model::point<T, 2, cs_type> point_type;
 173     typedef model::point<T, 3, cs::cartesian> vector_type;
 174
 175     collected_vector()
 176     {}
 177
 178     template <typename Point>
 179     collected_vector(Point const& p1, Point const& p2)
 180         : origin(get<0>(p1), get<1>(p1))
 181     {
 182         origin = detail::return_normalized<point_type>(origin);
 183
 184         using namespace geometry::formula;
 185         prev = sph_to_cart3d<vector_type>(p1);
 186         next = sph_to_cart3d<vector_type>(p2);
 187         direction = cross_product(prev, next);
 188     }
 189
 190     bool normalize()
 191     {
 192         T magnitude_sqr = dot_product(direction, direction);
 193
 194         // NOTE: shouldn't here math::equals() be called?
 195         if (magnitude_sqr > 0)
 196         {
 197             divide_value(direction, math::sqrt(magnitude_sqr));
 198             return true;
 199         }
 200
 201         return false;
 202     }
 203
 204     bool operator<(collected_vector const& other) const
 205     {
 206         if (math::equals(get<0>(origin), get<0>(other.origin)))
 207         {
 208             if (math::equals(get<1>(origin), get<1>(other.origin)))
 209             {
 210                 if (math::equals(get<0>(direction), get<0>(other.direction)))
 211                 {
 212                     if (math::equals(get<1>(direction), get<1>(other.direction)))
 213                     {
 214                         return get<2>(direction) < get<2>(other.direction);
 215                     }
 216                     return get<1>(direction) < get<1>(other.direction);
 217                 }
 218                 return get<0>(direction) < get<0>(other.direction);
 219             }
 220             return get<1>(origin) < get<1>(other.origin);
 221         }
 222         return get<0>(origin) < get<0>(other.origin);
 223     }
 224
 225     // For consistency with side and intersection strategies used by relops
 226     // IMPORTANT: this method should be called for previous vector
 227     // and next vector should be passed as parameter
 228     bool next_is_collinear(collected_vector const& other) const
 229     {
 230         return formula::sph_side_value(direction, other.next) == 0;
 231     }
 232
 233     // For std::equals
 234     bool operator==(collected_vector const& other) const
 235     {
 236         return math::equals(get<0>(origin), get<0>(other.origin))
 237             && math::equals(get<1>(origin), get<1>(other.origin))
 238             && is_collinear(other);
 239     }
 240
 241 private:
 242     // For consistency with side and intersection strategies used by relops
 243     bool is_collinear(collected_vector const& other) const
 244     {
 245         return formula::sph_side_value(direction, other.prev) == 0
 246             && formula::sph_side_value(direction, other.next) == 0;
 247     }
 248
 249     /*bool same_direction(collected_vector const& other) const
 250     {
 251         return math::equals_with_epsilon(get<0>(direction), get<0>(other.direction))
 252             && math::equals_with_epsilon(get<1>(direction), get<1>(other.direction))
 253             && math::equals_with_epsilon(get<2>(direction), get<2>(other.direction));
 254     }*/
 255
 256     point_type origin; // used for sorting and equality check
 257     vector_type direction; // used for sorting, only in operator<
 258     vector_type prev; // used for collinearity check, only in operator==
 259     vector_type next; // used for collinearity check
 260 };
 261
 262 // Specialization for spherical polar
 263 template <typename T, typename Geometry, typename CT>
 264 struct collected_vector
 265     <
 266         T, Geometry,
 267         strategy::side::spherical_side_formula<CT>,
 268         spherical_polar_tag
 269     >
 270     : public collected_vector
 271         <
 272             T, Geometry,
 273             strategy::side::spherical_side_formula<CT>,
 274             spherical_equatorial_tag
 275         >
 276 {
 277     typedef collected_vector
 278         <
 279             T, Geometry,
 280             strategy::side::spherical_side_formula<CT>,
 281             spherical_equatorial_tag
 282         > base_type;
 283
 284     collected_vector() {}
 285
 286     template <typename Point>
 287     collected_vector(Point const& p1, Point const& p2)
 288         : base_type(to_equatorial(p1), to_equatorial(p2))
 289     {}
 290
 291 private:
 292     template <typename Point>
 293     Point polar_to_equatorial(Point const& p)
 294     {
 295         typedef typename coordinate_type<Point>::type coord_type;
 296
 297         typedef math::detail::constants_on_spheroid
 298             <
 299                 coord_type,
 300                 typename coordinate_system<Point>::type::units
 301             > constants;
 302
 303         coord_type const pi_2 = constants::half_period() / 2;
 304
 305         Point res = p;
 306         set<1>(res, pi_2 - get<1>(p));
 307         return res;
 308     }
 309 };
 310
 311
 312 #ifndef DOXYGEN_NO_DETAIL
 313 namespace detail { namespace collect_vectors
 314 {
 315
 316
 317 template <typename Range, typename Collection>
 318 struct range_collect_vectors
 319 {
 320     typedef typename boost::range_value<Collection>::type item_type;
 321     typedef typename item_type::type calculation_type;
 322
 323     static inline void apply(Collection& collection, Range const& range)
 324     {
 325         typedef geometry::detail::normalized_view
 326             <
 327                 Range const
 328             > normalized_range_type;
 329
 330         apply_impl(collection, normalized_range_type(range));
 331     }
 332
 333 private:
 334     template <typename NormalizedRange>
 335     static inline void apply_impl(Collection& collection, NormalizedRange const& range)
 336     {
 337         if (boost::size(range) < 2)
 338         {
 339             return;
 340         }
 341
 342         typedef typename boost::range_size<Collection>::type collection_size_t;
 343         collection_size_t c_old_size = boost::size(collection);
 344
 345         typedef typename boost::range_iterator<NormalizedRange const>::type iterator;
 346
 347         bool is_first = true;
 348         iterator it = boost::begin(range);
 349
 350         for (iterator prev = it++;
 351             it != boost::end(range);
 352             prev = it++)
 353         {
 354             typename boost::range_value<Collection>::type v(*prev, *it);
 355
 356             // Normalize the vector -> this results in points+direction
 357             // and is comparible between geometries
 358             // Avoid non-duplicate points (AND division by zero)
 359             if (v.normalize())
 360             {
 361                 // Avoid non-direction changing points
 362                 if (is_first || ! collection.back().next_is_collinear(v))
 363                 {
 364                     collection.push_back(v);
 365                 }
 366                 is_first = false;
 367             }
 368         }
 369
 370         // If first one has same direction as last one, remove first one
 371         collection_size_t collected_count = boost::size(collection) - c_old_size;
 372         if ( collected_count > 1 )
 373         {
 374             typedef typename boost::range_iterator<Collection>::type c_iterator;
 375             c_iterator first = range::pos(collection, c_old_size);
 376
 377             if (collection.back().next_is_collinear(*first) )
 378             {
 379                 //collection.erase(first);
 380                 // O(1) instead of O(N)
 381                 *first = collection.back();
 382                 collection.pop_back();
 383             }
 384         }
 385     }
 386 };
 387
 388
 389 // Default version (cartesian)
 390 template <typename Box, typename Collection, typename CSTag = typename cs_tag<Box>::type>
 391 struct box_collect_vectors
 392 {
 393     // Calculate on coordinate type, but if it is integer,
 394     // then use double
 395     typedef typename boost::range_value<Collection>::type item_type;
 396     typedef typename item_type::type calculation_type;
 397
 398     static inline void apply(Collection& collection, Box const& box)
 399     {
 400         typename point_type<Box>::type lower_left, lower_right,
 401             upper_left, upper_right;
 402         geometry::detail::assign_box_corners(box, lower_left, lower_right,
 403             upper_left, upper_right);
 404
 405         typedef typename boost::range_value<Collection>::type item;
 406
 407         collection.push_back(item(get<0>(lower_left), get<1>(lower_left), 0, 1));
 408         collection.push_back(item(get<0>(upper_left), get<1>(upper_left), 1, 0));
 409         collection.push_back(item(get<0>(upper_right), get<1>(upper_right), 0, -1));
 410         collection.push_back(item(get<0>(lower_right), get<1>(lower_right), -1, 0));
 411     }
 412 };
 413
 414 // NOTE: This is not fully correct because Box in spherical and geographic
 415 // cordinate systems cannot be seen as Polygon
 416 template <typename Box, typename Collection>
 417 struct box_collect_vectors<Box, Collection, spherical_equatorial_tag>
 418 {
 419     static inline void apply(Collection& collection, Box const& box)
 420     {
 421         typename point_type<Box>::type lower_left, lower_right,
 422                 upper_left, upper_right;
 423         geometry::detail::assign_box_corners(box, lower_left, lower_right,
 424                 upper_left, upper_right);
 425
 426         typedef typename boost::range_value<Collection>::type item;
 427
 428         collection.push_back(item(lower_left, upper_left));
 429         collection.push_back(item(upper_left, upper_right));
 430         collection.push_back(item(upper_right, lower_right));
 431         collection.push_back(item(lower_right, lower_left));
 432     }
 433 };
 434
 435 template <typename Box, typename Collection>
 436 struct box_collect_vectors<Box, Collection, spherical_polar_tag>
 437     : box_collect_vectors<Box, Collection, spherical_equatorial_tag>
 438 {};
 439
 440 template <typename Box, typename Collection>
 441 struct box_collect_vectors<Box, Collection, geographic_tag>
 442     : box_collect_vectors<Box, Collection, spherical_equatorial_tag>
 443 {};
 444
 445
 446 template <typename Polygon, typename Collection>
 447 struct polygon_collect_vectors
 448 {
 449     static inline void apply(Collection& collection, Polygon const& polygon)
 450     {
 451         typedef typename geometry::ring_type<Polygon>::type ring_type;
 452
 453         typedef range_collect_vectors<ring_type, Collection> per_range;
 454         per_range::apply(collection, exterior_ring(polygon));
 455
 456         typename interior_return_type<Polygon const>::type
 457             rings = interior_rings(polygon);
 458         for (typename detail::interior_iterator<Polygon const>::type
 459                 it = boost::begin(rings); it != boost::end(rings); ++it)
 460         {
 461             per_range::apply(collection, *it);
 462         }
 463     }
 464 };
 465
 466
 467 template <typename MultiGeometry, typename Collection, typename SinglePolicy>
 468 struct multi_collect_vectors
 469 {
 470     static inline void apply(Collection& collection, MultiGeometry const& multi)
 471     {
 472         for (typename boost::range_iterator<MultiGeometry const>::type
 473                 it = boost::begin(multi);
 474             it != boost::end(multi);
 475             ++it)
 476         {
 477             SinglePolicy::apply(collection, *it);
 478         }
 479     }
 480 };
 481
 482
 483 }} // namespace detail::collect_vectors
 484 #endif // DOXYGEN_NO_DETAIL
 485
 486
 487
 488 #ifndef DOXYGEN_NO_DISPATCH
 489 namespace dispatch
 490 {
 491
 492
 493 template
 494 <
 495     typename Tag,
 496     typename Collection,
 497     typename Geometry
 498 >
 499 struct collect_vectors
 500 {
 501     static inline void apply(Collection&, Geometry const&)
 502     {}
 503 };
 504
 505
 506 template <typename Collection, typename Box>
 507 struct collect_vectors<box_tag, Collection, Box>
 508     : detail::collect_vectors::box_collect_vectors<Box, Collection>
 509 {};
 510
 511
 512
 513 template <typename Collection, typename Ring>
 514 struct collect_vectors<ring_tag, Collection, Ring>
 515     : detail::collect_vectors::range_collect_vectors<Ring, Collection>
 516 {};
 517
 518
 519 template <typename Collection, typename LineString>
 520 struct collect_vectors<linestring_tag, Collection, LineString>
 521     : detail::collect_vectors::range_collect_vectors<LineString, Collection>
 522 {};
 523
 524
 525 template <typename Collection, typename Polygon>
 526 struct collect_vectors<polygon_tag, Collection, Polygon>
 527     : detail::collect_vectors::polygon_collect_vectors<Polygon, Collection>
 528 {};
 529
 530
 531 template <typename Collection, typename MultiPolygon>
 532 struct collect_vectors<multi_polygon_tag, Collection, MultiPolygon>
 533     : detail::collect_vectors::multi_collect_vectors
 534         <
 535             MultiPolygon,
 536             Collection,
 537             detail::collect_vectors::polygon_collect_vectors
 538             <
 539                 typename boost::range_value<MultiPolygon>::type,
 540                 Collection
 541             >
 542         >
 543 {};
 544
 545
 546
 547 } // namespace dispatch
 548 #endif
 549
 550
 551 /*!
 552     \ingroup collect_vectors
 553     \tparam Collection Collection type, should be e.g. std::vector<>
 554     \tparam Geometry geometry type
 555     \param collection the collection of vectors
 556     \param geometry the geometry to make collect_vectors
 557 */
 558 template <typename Collection, typename Geometry>
 559 inline void collect_vectors(Collection& collection, Geometry const& geometry)
 560 {
 561     concepts::check<Geometry const>();
 562
 563     dispatch::collect_vectors
 564         <
 565             typename tag<Geometry>::type,
 566             Collection,
 567             Geometry
 568         >::apply(collection, geometry);
 569 }
 570
 571
 572 }} // namespace boost::geometry
 573
 574
 575 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EQUALS_COLLECT_VECTORS_HPP