geometry/algorithms/merge_elements.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2022-2023, Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Licensed under the Boost Software License version 1.0.
0009 // http://www.boost.org/users/license.html
0010
0011 #ifndef BOOST_GEOMETRY_ALGORITHMS_MERGE_ELEMENTS_HPP
0012 #define BOOST_GEOMETRY_ALGORITHMS_MERGE_ELEMENTS_HPP
0013
0014
0015 #include <vector>
0016
0017 #include <boost/geometry/algorithms/detail/point_on_border.hpp>
0018 #include <boost/geometry/algorithms/detail/visit.hpp>
0019 #include <boost/geometry/algorithms/difference.hpp>
0020 #include <boost/geometry/algorithms/union.hpp>
0021 #include <boost/geometry/core/coordinate_system.hpp>
0022 #include <boost/geometry/core/coordinate_type.hpp>
0023 #include <boost/geometry/core/point_type.hpp>
0024 #include <boost/geometry/core/tag.hpp>
0025 #include <boost/geometry/core/tags.hpp>
0026 #include <boost/geometry/core/visit.hpp>
0027 #include <boost/geometry/geometries/point.hpp>
0028 #include <boost/geometry/policies/compare.hpp>
0029 #include <boost/geometry/util/range.hpp>
0030 #include <boost/geometry/strategies/relate/cartesian.hpp>
0031 #include <boost/geometry/strategies/relate/geographic.hpp>
0032 #include <boost/geometry/strategies/relate/spherical.hpp>
0033
0034
0035 namespace boost { namespace geometry
0036 {
0037
0038
0039 #ifndef DOXYGEN_NO_DETAIL
0040 namespace detail { namespace merge_elements
0041 {
0042
0043
0044 template <typename T>
0045 using is_pla = util::bool_constant<util::is_pointlike<T>::value || util::is_linear<T>::value || util::is_areal<T>::value>;
0046 template <typename T, typename ...Ts>
0047 struct are_areal : util::bool_constant<util::is_areal<T>::value && are_areal<Ts...>::value> {};
0048 template <typename T>
0049 struct are_areal<T> : util::is_areal<T> {};
0050 template <typename T, typename ...Ts>
0051 struct are_linear : util::bool_constant<util::is_linear<T>::value && are_linear<Ts...>::value> {};
0052 template <typename T>
0053 struct are_linear<T> : util::is_linear<T> {};
0054 template <typename T, typename ...Ts>
0055 struct are_pointlike : util::bool_constant<util::is_pointlike<T>::value && are_pointlike<Ts...>::value> {};
0056 template <typename T>
0057 struct are_pointlike<T> : util::is_pointlike<T> {};
0058 template <typename ...Ts>
0059 using are_same_kind = util::bool_constant<are_areal<Ts...>::value || are_linear<Ts...>::value || are_pointlike<Ts...>::value>;
0060
0061
0062 template
0063 <
0064     typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
0065     std::enable_if_t<util::is_areal<Geometry>::value, int> = 0
0066 >
0067 inline void distribute_element(Geometry const& geometry, It it, PointLike& , Linear&, Areal& areal)
0068 {
0069     typename geometry::point_type<Geometry>::type point;
0070     if (geometry::point_on_border(point, geometry))
0071     {
0072         using point_t = typename Areal::value_type::first_type;
0073         areal.emplace_back(point_t(geometry::get<0>(point), geometry::get<1>(point)), it);
0074     }
0075 }
0076
0077 template
0078 <
0079     typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
0080     std::enable_if_t<util::is_linear<Geometry>::value, int> = 0
0081 >
0082 inline void distribute_element(Geometry const& geometry, It it, PointLike& , Linear& linear, Areal& )
0083 {
0084     typename geometry::point_type<Geometry>::type point;
0085     if (geometry::point_on_border(point, geometry))
0086     {
0087         using point_t = typename Linear::value_type::first_type;
0088         linear.emplace_back(point_t(geometry::get<0>(point), geometry::get<1>(point)), it);
0089     }
0090 }
0091
0092 template
0093 <
0094     typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
0095     std::enable_if_t<util::is_pointlike<Geometry>::value, int> = 0
0096 >
0097 inline void distribute_element(Geometry const& geometry, It it, PointLike& pointlike, Linear& , Areal& )
0098 {
0099     typename geometry::point_type<Geometry>::type point;
0100     if (geometry::point_on_border(point, geometry))
0101     {
0102         using point_t = typename Linear::value_type::first_type;
0103         pointlike.emplace_back(point_t(geometry::get<0>(point), geometry::get<1>(point)), it);
0104     }
0105 }
0106
0107 template
0108 <
0109     typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
0110     std::enable_if_t<! is_pla<Geometry>::value, int> = 0
0111 >
0112 inline void distribute_element(Geometry const& , It const&, PointLike const& , Linear const&, Areal const&)
0113 {}
0114
0115
0116 template
0117 <
0118     typename Geometry, typename MultiGeometry,
0119     std::enable_if_t<are_same_kind<Geometry, MultiGeometry>::value, int> = 0
0120 >
0121 inline void convert(Geometry const& geometry, MultiGeometry& result)
0122 {
0123     geometry::convert(geometry, result);
0124 }
0125
0126 template
0127 <
0128     typename Geometry, typename MultiGeometry,
0129     std::enable_if_t<! are_same_kind<Geometry, MultiGeometry>::value, int> = 0
0130 >
0131 inline void convert(Geometry const& , MultiGeometry const& )
0132 {}
0133
0134
0135 template
0136 <
0137     typename Geometry1, typename Geometry2, typename MultiGeometry, typename Strategy,
0138     std::enable_if_t<are_same_kind<Geometry1, Geometry2, MultiGeometry>::value, int> = 0
0139 >
0140 inline void union_(Geometry1 const& geometry1, Geometry2 const& geometry2, MultiGeometry& result, Strategy const& strategy)
0141 {
0142     geometry::union_(geometry1, geometry2, result, strategy);
0143 }
0144
0145 template
0146 <
0147     typename Geometry1, typename Geometry2, typename MultiGeometry, typename Strategy,
0148     std::enable_if_t<! are_same_kind<Geometry1, Geometry2, MultiGeometry>::value, int> = 0
0149 >
0150 inline void union_(Geometry1 const& , Geometry2 const& , MultiGeometry const& , Strategy const&)
0151 {}
0152
0153
0154 template <typename It>
0155 struct merge_data
0156 {
0157     merge_data(It first_, It last_)
0158         : first(first_), last(last_)
0159     {}
0160     It first, last;
0161     bool merge_results = false;
0162 };
0163
0164 template <typename GeometryCollection, typename RandomIt, typename MultiGeometry, typename Strategy>
0165 inline void merge(RandomIt const first, RandomIt const last, MultiGeometry& out, Strategy const& strategy)
0166 {
0167     auto const size = last - first;
0168     if (size <= 0)
0169     {
0170         return;
0171     }
0172
0173     auto const less = [](auto const& l, auto const& r)
0174     {
0175         return geometry::less<void, -1, Strategy>()(l.first, r.first);
0176     };
0177
0178     std::vector<merge_data<RandomIt>> stack_in;
0179     std::vector<MultiGeometry> stack_out;
0180     stack_in.reserve(size / 2 + 1);
0181     stack_out.reserve(size / 2 + 1);
0182
0183     stack_in.emplace_back(first, last);
0184
0185     while (! stack_in.empty())
0186     {
0187         auto & b = stack_in.back();
0188         if (! b.merge_results)
0189         {
0190             auto const s = b.last - b.first;
0191             if (s > 2)
0192             {
0193                 RandomIt const mid = b.first + s / 2;
0194                 std::nth_element(b.first, mid, b.last, less);
0195                 RandomIt const fir = b.first;
0196                 RandomIt const las = b.last;
0197                 b.merge_results = true;
0198                 stack_in.emplace_back(fir, mid);
0199                 stack_in.emplace_back(mid, las);
0200             }
0201             else if (s == 2)
0202             {
0203                 MultiGeometry result;
0204                 // VERSION 1
0205 //                traits::iter_visit<GeometryCollection>::apply([&](auto const& g1)
0206 //                {
0207 //                    traits::iter_visit<GeometryCollection>::apply([&](auto const& g2)
0208 //                    {
0209 //                        merge_elements::union_(g1, g2, result, strategy);
0210 //                    }, (b.first + 1)->second);
0211 //                }, b.first->second);
0212                 // VERSION 2
0213                 // calling iter_visit non-recursively seems to decrease compilation time
0214                 // greately with GCC
0215                 MultiGeometry temp1, temp2;
0216                 traits::iter_visit<GeometryCollection>::apply([&](auto const& g1)
0217                 {
0218                     merge_elements::convert(g1, temp1);
0219                 }, b.first->second);
0220                 traits::iter_visit<GeometryCollection>::apply([&](auto const& g2)
0221                 {
0222                     merge_elements::convert(g2, temp2);
0223                 }, (b.first + 1)->second);
0224                 geometry::union_(temp1, temp2, result, strategy);
0225
0226                 stack_out.push_back(std::move(result));
0227                 stack_in.pop_back();
0228             }
0229             else if (s == 1)
0230             {
0231                 MultiGeometry result;
0232                 traits::iter_visit<GeometryCollection>::apply([&](auto const& g)
0233                 {
0234                     merge_elements::convert(g, result);
0235                 }, b.first->second);
0236                 stack_out.push_back(std::move(result));
0237                 stack_in.pop_back();
0238             }
0239         }
0240         else if (b.merge_results)
0241         {
0242             MultiGeometry m2 = std::move(stack_out.back());
0243             stack_out.pop_back();
0244             MultiGeometry m1 = std::move(stack_out.back());
0245             stack_out.pop_back();
0246             MultiGeometry result;
0247             geometry::union_(m1, m2, result, strategy);
0248             stack_out.push_back(std::move(result));
0249             stack_in.pop_back();
0250         }
0251     }
0252
0253     out = std::move(stack_out.back());
0254 }
0255
0256 template <typename MultiGeometry, typename Geometry, typename Strategy>
0257 inline void subtract(MultiGeometry & multi, Geometry const& geometry, Strategy const& strategy)
0258 {
0259     MultiGeometry temp;
0260     geometry::difference(multi, geometry, temp, strategy);
0261     multi = std::move(temp);
0262 }
0263
0264 struct merge_gc
0265 {
0266     template <typename GeometryCollection, typename Strategy>
0267     static void apply(GeometryCollection const& geometry_collection,
0268                       GeometryCollection & out,
0269                       Strategy const& strategy)
0270     {
0271         using original_point_t = typename geometry::point_type<GeometryCollection>::type;
0272         using iterator_t = typename boost::range_iterator<GeometryCollection const>::type;
0273         using coordinate_t = typename geometry::coordinate_type<original_point_t>::type;
0274         using cs_t = typename geometry::coordinate_system<original_point_t>::type;
0275         using point_t = model::point<coordinate_t, 2, cs_t>;
0276
0277         using multi_point_t = typename util::sequence_find_if
0278             <
0279                 typename traits::geometry_types<std::remove_const_t<GeometryCollection>>::type,
0280                 util::is_multi_point
0281             >::type;
0282         using multi_linestring_t = typename util::sequence_find_if
0283             <
0284                 typename traits::geometry_types<std::remove_const_t<GeometryCollection>>::type,
0285                 util::is_multi_linestring
0286             >::type;
0287         using multi_polygon_t = typename util::sequence_find_if
0288             <
0289                 typename traits::geometry_types<std::remove_const_t<GeometryCollection>>::type,
0290                 util::is_multi_polygon
0291             >::type;
0292
0293         // NOTE: Right now GC containing all of the above is required but technically
0294         //       we could allow only some combinations and the algorithm below could
0295         //       normalize GC accordingly.
0296
0297         multi_point_t multi_point;
0298         multi_linestring_t multi_linestring;
0299         multi_polygon_t multi_polygon;
0300
0301         std::vector<std::pair<point_t, iterator_t>> pointlike;
0302         std::vector<std::pair<point_t, iterator_t>> linear;
0303         std::vector<std::pair<point_t, iterator_t>> areal;
0304
0305         detail::visit_breadth_first_impl<true>::apply([&](auto const& g, auto it)
0306         {
0307             merge_elements::distribute_element(g, it, pointlike, linear, areal);
0308             return true;
0309         }, geometry_collection);
0310
0311         // TODO: make this optional?
0312         // TODO: merge linear at the end? (difference can break linear rings, their parts would be joined)
0313         merge<GeometryCollection>(pointlike.begin(), pointlike.end(), multi_point, strategy);
0314         merge<GeometryCollection>(linear.begin(), linear.end(), multi_linestring, strategy);
0315         merge<GeometryCollection>(areal.begin(), areal.end(), multi_polygon, strategy);
0316
0317         // L \ A
0318         subtract(multi_linestring, multi_polygon, strategy);
0319         // P \ A
0320         subtract(multi_point, multi_polygon, strategy);
0321         // P \ L
0322         subtract(multi_point, multi_linestring, strategy);
0323
0324         if (! geometry::is_empty(multi_point))
0325         {
0326             range::emplace_back(out, std::move(multi_point));
0327         }
0328         if (! geometry::is_empty(multi_linestring))
0329         {
0330             range::emplace_back(out, std::move(multi_linestring));
0331         }
0332         if (! geometry::is_empty(multi_polygon))
0333         {
0334             range::emplace_back(out, std::move(multi_polygon));
0335         }
0336     }
0337 };
0338
0339
0340 }} // namespace detail::merge_elements
0341 #endif // DOXYGEN_NO_DETAIL
0342
0343
0344 #ifndef DOXYGEN_NO_DISPATCH
0345 namespace dispatch
0346 {
0347
0348
0349 template <typename Geometry, typename Tag = typename tag<Geometry>::type>
0350 struct merge_elements
0351     : not_implemented<Geometry, Tag>
0352 {};
0353
0354 template <typename GeometryCollection>
0355 struct merge_elements<GeometryCollection, geometry_collection_tag>
0356     : geometry::detail::merge_elements::merge_gc
0357 {};
0358
0359
0360 } // namespace dispatch
0361 #endif
0362
0363
0364 namespace resolve_strategy
0365 {
0366
0367 template <typename Strategy>
0368 struct merge_elements
0369 {
0370     template <typename Geometry>
0371     static void apply(Geometry const& geometry, Geometry & out, Strategy const& strategy)
0372     {
0373         dispatch::merge_elements
0374             <
0375                 Geometry
0376             >::apply(geometry, out, strategy);
0377     }
0378 };
0379
0380 template <>
0381 struct merge_elements<default_strategy>
0382 {
0383     template <typename Geometry>
0384     static void apply(Geometry const& geometry, Geometry & out, default_strategy)
0385     {
0386         using strategy_type = typename strategies::relate::services::default_strategy
0387             <
0388                 Geometry, Geometry
0389             >::type;
0390
0391         dispatch::merge_elements
0392             <
0393                 Geometry
0394             >::apply(geometry, out, strategy_type());
0395     }
0396 };
0397
0398 } // namespace resolve_strategy
0399
0400
0401 template <typename Geometry, typename Strategy>
0402 inline void merge_elements(Geometry const& geometry, Geometry & out, Strategy const& strategy)
0403 {
0404     resolve_strategy::merge_elements
0405         <
0406             Strategy
0407         >::apply(geometry, out, strategy);
0408 }
0409
0410
0411 template <typename Geometry>
0412 inline void merge_elements(Geometry const& geometry, Geometry & out)
0413 {
0414     resolve_strategy::merge_elements
0415         <
0416             default_strategy
0417         >::apply(geometry, out, default_strategy());
0418 }
0419
0420
0421 }} // namespace boost::geometry
0422
0423
0424 #endif // BOOST_GEOMETRY_ALGORITHMS_MERGE_ELEMENTS_HPP