425 lines
14 KiB
C++
425 lines
14 KiB
C++
|
// Boost.Geometry
|
||
|
|
|
||
|
|
// Copyright (c) 2022, Oracle and/or its affiliates.
|
||
|
|
|
||
|
|
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
|
||
|
|
|
||
|
|
// Licensed under the Boost Software License version 1.0.
|
||
|
|
// http://www.boost.org/users/license.html
|
||
|
|
|
||
|
|
#ifndef BOOST_GEOMETRY_ALGORITHMS_MERGE_ELEMENTS_HPP
|
||
|
|
#define BOOST_GEOMETRY_ALGORITHMS_MERGE_ELEMENTS_HPP
|
||
|
|
|
||
|
|
|
||
|
|
#include <vector>
|
||
|
|
|
||
|
|
#include <boost/geometry/algorithms/detail/point_on_border.hpp>
|
||
|
|
#include <boost/geometry/algorithms/detail/visit.hpp>
|
||
|
|
#include <boost/geometry/algorithms/difference.hpp>
|
||
|
|
#include <boost/geometry/algorithms/union.hpp>
|
||
|
|
#include <boost/geometry/core/coordinate_system.hpp>
|
||
|
|
#include <boost/geometry/core/coordinate_type.hpp>
|
||
|
|
#include <boost/geometry/core/point_type.hpp>
|
||
|
|
#include <boost/geometry/core/tag.hpp>
|
||
|
|
#include <boost/geometry/core/tags.hpp>
|
||
|
|
#include <boost/geometry/core/visit.hpp>
|
||
|
|
#include <boost/geometry/geometries/point.hpp>
|
||
|
|
#include <boost/geometry/policies/compare.hpp>
|
||
|
|
#include <boost/geometry/util/range.hpp>
|
||
|
|
#include <boost/geometry/strategies/relate/cartesian.hpp>
|
||
|
|
#include <boost/geometry/strategies/relate/geographic.hpp>
|
||
|
|
#include <boost/geometry/strategies/relate/spherical.hpp>
|
||
|
|
#include <boost/geometry/strategies/spherical/compare.hpp>
|
||
|
|
|
||
|
|
|
||
|
|
namespace boost { namespace geometry
|
||
|
|
{
|
||
|
|
|
||
|
|
|
||
|
|
#ifndef DOXYGEN_NO_DETAIL
|
||
|
|
namespace detail { namespace merge_elements
|
||
|
|
{
|
||
|
|
|
||
|
|
|
||
|
|
template <typename T>
|
||
|
|
using is_pla = util::bool_constant<util::is_pointlike<T>::value || util::is_linear<T>::value || util::is_areal<T>::value>;
|
||
|
|
template <typename T, typename ...Ts>
|
||
|
|
struct are_areal : util::bool_constant<util::is_areal<T>::value && are_areal<Ts...>::value> {};
|
||
|
|
template <typename T>
|
||
|
|
struct are_areal<T> : util::is_areal<T> {};
|
||
|
|
template <typename T, typename ...Ts>
|
||
|
|
struct are_linear : util::bool_constant<util::is_linear<T>::value && are_linear<Ts...>::value> {};
|
||
|
|
template <typename T>
|
||
|
|
struct are_linear<T> : util::is_linear<T> {};
|
||
|
|
template <typename T, typename ...Ts>
|
||
|
|
struct are_pointlike : util::bool_constant<util::is_pointlike<T>::value && are_pointlike<Ts...>::value> {};
|
||
|
|
template <typename T>
|
||
|
|
struct are_pointlike<T> : util::is_pointlike<T> {};
|
||
|
|
template <typename ...Ts>
|
||
|
|
using are_same_kind = util::bool_constant<are_areal<Ts...>::value || are_linear<Ts...>::value || are_pointlike<Ts...>::value>;
|
||
|
|
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
|
||
|
|
std::enable_if_t<util::is_areal<Geometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void distribute_element(Geometry const& geometry, It it, PointLike& , Linear&, Areal& areal)
|
||
|
|
{
|
||
|
|
typename geometry::point_type<Geometry>::type point;
|
||
|
|
if (geometry::point_on_border(point, geometry))
|
||
|
|
{
|
||
|
|
using point_t = typename Areal::value_type::first_type;
|
||
|
|
areal.emplace_back(point_t(geometry::get<0>(point), geometry::get<1>(point)), it);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
|
||
|
|
std::enable_if_t<util::is_linear<Geometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void distribute_element(Geometry const& geometry, It it, PointLike& , Linear& linear, Areal& )
|
||
|
|
{
|
||
|
|
typename geometry::point_type<Geometry>::type point;
|
||
|
|
if (geometry::point_on_border(point, geometry))
|
||
|
|
{
|
||
|
|
using point_t = typename Linear::value_type::first_type;
|
||
|
|
linear.emplace_back(point_t(geometry::get<0>(point), geometry::get<1>(point)), it);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
|
||
|
|
std::enable_if_t<util::is_pointlike<Geometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void distribute_element(Geometry const& geometry, It it, PointLike& pointlike, Linear& , Areal& )
|
||
|
|
{
|
||
|
|
typename geometry::point_type<Geometry>::type point;
|
||
|
|
if (geometry::point_on_border(point, geometry))
|
||
|
|
{
|
||
|
|
using point_t = typename Linear::value_type::first_type;
|
||
|
|
pointlike.emplace_back(point_t(geometry::get<0>(point), geometry::get<1>(point)), it);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry, typename It, typename PointLike, typename Linear, typename Areal,
|
||
|
|
std::enable_if_t<! is_pla<Geometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void distribute_element(Geometry const& , It const&, PointLike const& , Linear const&, Areal const&)
|
||
|
|
{}
|
||
|
|
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry, typename MultiGeometry,
|
||
|
|
std::enable_if_t<are_same_kind<Geometry, MultiGeometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void convert(Geometry const& geometry, MultiGeometry& result)
|
||
|
|
{
|
||
|
|
geometry::convert(geometry, result);
|
||
|
|
}
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry, typename MultiGeometry,
|
||
|
|
std::enable_if_t<! are_same_kind<Geometry, MultiGeometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void convert(Geometry const& , MultiGeometry const& )
|
||
|
|
{}
|
||
|
|
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry1, typename Geometry2, typename MultiGeometry, typename Strategy,
|
||
|
|
std::enable_if_t<are_same_kind<Geometry1, Geometry2, MultiGeometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void union_(Geometry1 const& geometry1, Geometry2 const& geometry2, MultiGeometry& result, Strategy const& strategy)
|
||
|
|
{
|
||
|
|
geometry::union_(geometry1, geometry2, result, strategy);
|
||
|
|
}
|
||
|
|
|
||
|
|
template
|
||
|
|
<
|
||
|
|
typename Geometry1, typename Geometry2, typename MultiGeometry, typename Strategy,
|
||
|
|
std::enable_if_t<! are_same_kind<Geometry1, Geometry2, MultiGeometry>::value, int> = 0
|
||
|
|
>
|
||
|
|
inline void union_(Geometry1 const& , Geometry2 const& , MultiGeometry const& , Strategy const&)
|
||
|
|
{}
|
||
|
|
|
||
|
|
|
||
|
|
template <typename It>
|
||
|
|
struct merge_data
|
||
|
|
{
|
||
|
|
merge_data(It first_, It last_)
|
||
|
|
: first(first_), last(last_)
|
||
|
|
{}
|
||
|
|
It first, last;
|
||
|
|
bool merge_results = false;
|
||
|
|
};
|
||
|
|
|
||
|
|
template <typename GeometryCollection, typename RandomIt, typename MultiGeometry, typename Strategy>
|
||
|
|
inline void merge(RandomIt const first, RandomIt const last, MultiGeometry& out, Strategy const& strategy)
|
||
|
|
{
|
||
|
|
auto const size = last - first;
|
||
|
|
if (size <= 0)
|
||
|
|
{
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
auto const less = [](auto const& l, auto const& r)
|
||
|
|
{
|
||
|
|
return geometry::less<void, -1, typename Strategy::cs_tag>()(l.first, r.first);
|
||
|
|
};
|
||
|
|
|
||
|
|
std::vector<merge_data<RandomIt>> stack_in;
|
||
|
|
std::vector<MultiGeometry> stack_out;
|
||
|
|
stack_in.reserve(size / 2 + 1);
|
||
|
|
stack_out.reserve(size / 2 + 1);
|
||
|
|
|
||
|
|
stack_in.emplace_back(first, last);
|
||
|
|
|
||
|
|
while (! stack_in.empty())
|
||
|
|
{
|
||
|
|
auto & b = stack_in.back();
|
||
|
|
if (! b.merge_results)
|
||
|
|
{
|
||
|
|
auto const s = b.last - b.first;
|
||
|
|
if (s > 2)
|
||
|
|
{
|
||
|
|
RandomIt const mid = b.first + s / 2;
|
||
|
|
std::nth_element(b.first, mid, b.last, less);
|
||
|
|
RandomIt const fir = b.first;
|
||
|
|
RandomIt const las = b.last;
|
||
|
|
b.merge_results = true;
|
||
|
|
stack_in.emplace_back(fir, mid);
|
||
|
|
stack_in.emplace_back(mid, las);
|
||
|
|
}
|
||
|
|
else if (s == 2)
|
||
|
|
{
|
||
|
|
MultiGeometry result;
|
||
|
|
// VERSION 1
|
||
|
|
// traits::iter_visit<GeometryCollection>::apply([&](auto const& g1)
|
||
|
|
// {
|
||
|
|
// traits::iter_visit<GeometryCollection>::apply([&](auto const& g2)
|
||
|
|
// {
|
||
|
|
// merge_elements::union_(g1, g2, result, strategy);
|
||
|
|
// }, (b.first + 1)->second);
|
||
|
|
// }, b.first->second);
|
||
|
|
// VERSION 2
|
||
|
|
// calling iter_visit non-recursively seems to decrease compilation time
|
||
|
|
// greately with GCC
|
||
|
|
MultiGeometry temp1, temp2;
|
||
|
|
traits::iter_visit<GeometryCollection>::apply([&](auto const& g1)
|
||
|
|
{
|
||
|
|
merge_elements::convert(g1, temp1);
|
||
|
|
}, b.first->second);
|
||
|
|
traits::iter_visit<GeometryCollection>::apply([&](auto const& g2)
|
||
|
|
{
|
||
|
|
merge_elements::convert(g2, temp2);
|
||
|
|
}, (b.first + 1)->second);
|
||
|
|
geometry::union_(temp1, temp2, result, strategy);
|
||
|
|
|
||
|
|
stack_out.push_back(std::move(result));
|
||
|
|
stack_in.pop_back();
|
||
|
|
}
|
||
|
|
else if (s == 1)
|
||
|
|
{
|
||
|
|
MultiGeometry result;
|
||
|
|
traits::iter_visit<GeometryCollection>::apply([&](auto const& g)
|
||
|
|
{
|
||
|
|
merge_elements::convert(g, result);
|
||
|
|
}, b.first->second);
|
||
|
|
stack_out.push_back(std::move(result));
|
||
|
|
stack_in.pop_back();
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else if (b.merge_results)
|
||
|
|
{
|
||
|
|
MultiGeometry m2 = std::move(stack_out.back());
|
||
|
|
stack_out.pop_back();
|
||
|
|
MultiGeometry m1 = std::move(stack_out.back());
|
||
|
|
stack_out.pop_back();
|
||
|
|
MultiGeometry result;
|
||
|
|
geometry::union_(m1, m2, result, strategy);
|
||
|
|
stack_out.push_back(std::move(result));
|
||
|
|
stack_in.pop_back();
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
out = std::move(stack_out.back());
|
||
|
|
}
|
||
|
|
|
||
|
|
template <typename MultiGeometry, typename Geometry, typename Strategy>
|
||
|
|
inline void subtract(MultiGeometry & multi, Geometry const& geometry, Strategy const& strategy)
|
||
|
|
{
|
||
|
|
MultiGeometry temp;
|
||
|
|
geometry::difference(multi, geometry, temp, strategy);
|
||
|
|
multi = std::move(temp);
|
||
|
|
}
|
||
|
|
|
||
|
|
struct merge_gc
|
||
|
|
{
|
||
|
|
template <typename GeometryCollection, typename Strategy>
|
||
|
|
static void apply(GeometryCollection const& geometry_collection,
|
||
|
|
GeometryCollection & out,
|
||
|
|
Strategy const& strategy)
|
||
|
|
{
|
||
|
|
using original_point_t = typename geometry::point_type<GeometryCollection>::type;
|
||
|
|
using iterator_t = typename boost::range_iterator<GeometryCollection const>::type;
|
||
|
|
using coordinate_t = typename geometry::coordinate_type<original_point_t>::type;
|
||
|
|
using cs_t = typename geometry::coordinate_system<original_point_t>::type;
|
||
|
|
using point_t = model::point<coordinate_t, 2, cs_t>;
|
||
|
|
|
||
|
|
using multi_point_t = typename util::sequence_find_if
|
||
|
|
<
|
||
|
|
typename traits::geometry_types<std::remove_const_t<GeometryCollection>>::type,
|
||
|
|
util::is_multi_point
|
||
|
|
>::type;
|
||
|
|
using multi_linestring_t = typename util::sequence_find_if
|
||
|
|
<
|
||
|
|
typename traits::geometry_types<std::remove_const_t<GeometryCollection>>::type,
|
||
|
|
util::is_multi_linestring
|
||
|
|
>::type;
|
||
|
|
using multi_polygon_t = typename util::sequence_find_if
|
||
|
|
<
|
||
|
|
typename traits::geometry_types<std::remove_const_t<GeometryCollection>>::type,
|
||
|
|
util::is_multi_polygon
|
||
|
|
>::type;
|
||
|
|
|
||
|
|
// NOTE: Right now GC containing all of the above is required but technically
|
||
|
|
// we could allow only some combinations and the algorithm below could
|
||
|
|
// normalize GC accordingly.
|
||
|
|
|
||
|
|
multi_point_t multi_point;
|
||
|
|
multi_linestring_t multi_linestring;
|
||
|
|
multi_polygon_t multi_polygon;
|
||
|
|
|
||
|
|
std::vector<std::pair<point_t, iterator_t>> pointlike;
|
||
|
|
std::vector<std::pair<point_t, iterator_t>> linear;
|
||
|
|
std::vector<std::pair<point_t, iterator_t>> areal;
|
||
|
|
|
||
|
|
detail::visit_breadth_first_impl<true>::apply([&](auto const& g, auto it)
|
||
|
|
{
|
||
|
|
merge_elements::distribute_element(g, it, pointlike, linear, areal);
|
||
|
|
return true;
|
||
|
|
}, geometry_collection);
|
||
|
|
|
||
|
|
// TODO: make this optional?
|
||
|
|
// TODO: merge linear at the end? (difference can break linear rings, their parts would be joined)
|
||
|
|
merge<GeometryCollection>(pointlike.begin(), pointlike.end(), multi_point, strategy);
|
||
|
|
merge<GeometryCollection>(linear.begin(), linear.end(), multi_linestring, strategy);
|
||
|
|
merge<GeometryCollection>(areal.begin(), areal.end(), multi_polygon, strategy);
|
||
|
|
|
||
|
|
// L \ A
|
||
|
|
subtract(multi_linestring, multi_polygon, strategy);
|
||
|
|
// P \ A
|
||
|
|
subtract(multi_point, multi_polygon, strategy);
|
||
|
|
// P \ L
|
||
|
|
subtract(multi_point, multi_linestring, strategy);
|
||
|
|
|
||
|
|
if (! geometry::is_empty(multi_point))
|
||
|
|
{
|
||
|
|
range::emplace_back(out, std::move(multi_point));
|
||
|
|
}
|
||
|
|
if (! geometry::is_empty(multi_linestring))
|
||
|
|
{
|
||
|
|
range::emplace_back(out, std::move(multi_linestring));
|
||
|
|
}
|
||
|
|
if (! geometry::is_empty(multi_polygon))
|
||
|
|
{
|
||
|
|
range::emplace_back(out, std::move(multi_polygon));
|
||
|
|
}
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
|
||
|
|
}} // namespace detail::merge_elements
|
||
|
|
#endif // DOXYGEN_NO_DETAIL
|
||
|
|
|
||
|
|
|
||
|
|
#ifndef DOXYGEN_NO_DISPATCH
|
||
|
|
namespace dispatch
|
||
|
|
{
|
||
|
|
|
||
|
|
|
||
|
|
template <typename Geometry, typename Tag = typename tag<Geometry>::type>
|
||
|
|
struct merge_elements
|
||
|
|
: not_implemented<Geometry, Tag>
|
||
|
|
{};
|
||
|
|
|
||
|
|
template <typename GeometryCollection>
|
||
|
|
struct merge_elements<GeometryCollection, geometry_collection_tag>
|
||
|
|
: geometry::detail::merge_elements::merge_gc
|
||
|
|
{};
|
||
|
|
|
||
|
|
|
||
|
|
} // namespace dispatch
|
||
|
|
#endif
|
||
|
|
|
||
|
|
|
||
|
|
namespace resolve_strategy
|
||
|
|
{
|
||
|
|
|
||
|
|
template <typename Strategy>
|
||
|
|
struct merge_elements
|
||
|
|
{
|
||
|
|
template <typename Geometry>
|
||
|
|
static void apply(Geometry const& geometry, Geometry & out, Strategy const& strategy)
|
||
|
|
{
|
||
|
|
dispatch::merge_elements
|
||
|
|
<
|
||
|
|
Geometry
|
||
|
|
>::apply(geometry, out, strategy);
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
template <>
|
||
|
|
struct merge_elements<default_strategy>
|
||
|
|
{
|
||
|
|
template <typename Geometry>
|
||
|
|
static void apply(Geometry const& geometry, Geometry & out, default_strategy)
|
||
|
|
{
|
||
|
|
using strategy_type = typename strategies::relate::services::default_strategy
|
||
|
|
<
|
||
|
|
Geometry, Geometry
|
||
|
|
>::type;
|
||
|
|
|
||
|
|
dispatch::merge_elements
|
||
|
|
<
|
||
|
|
Geometry
|
||
|
|
>::apply(geometry, out, strategy_type());
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
} // namespace resolve_strategy
|
||
|
|
|
||
|
|
|
||
|
|
template <typename Geometry, typename Strategy>
|
||
|
|
inline void merge_elements(Geometry const& geometry, Geometry & out, Strategy const& strategy)
|
||
|
|
{
|
||
|
|
resolve_strategy::merge_elements
|
||
|
|
<
|
||
|
|
Strategy
|
||
|
|
>::apply(geometry, out, strategy);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
template <typename Geometry>
|
||
|
|
inline void merge_elements(Geometry const& geometry, Geometry & out)
|
||
|
|
{
|
||
|
|
resolve_strategy::merge_elements
|
||
|
|
<
|
||
|
|
default_strategy
|
||
|
|
>::apply(geometry, out, default_strategy());
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
}} // namespace boost::geometry
|
||
|
|
|
||
|
|
|
||
|
|
#endif // BOOST_GEOMETRY_ALGORITHMS_MERGE_ELEMENTS_HPP
|