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