From 1df1e7021e16735b6d471a77b78d63777f34b435 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Wed, 19 Oct 2022 11:43:48 +0200 Subject: [PATCH] [buffer] harmonize geographic strategy code and share code --- .../geographic/buffer_end_round.hpp | 42 +++----- .../strategies/geographic/buffer_helper.hpp | 101 ++++++++++++++++++ .../geographic/buffer_join_miter.hpp | 54 ++-------- .../geographic/buffer_join_round.hpp | 56 +++------- .../geographic/buffer_point_circle.hpp | 66 ++++++------ .../geographic/buffer_side_straight.hpp | 49 +++------ 6 files changed, 179 insertions(+), 189 deletions(-) create mode 100644 include/boost/geometry/strategies/geographic/buffer_helper.hpp diff --git a/include/boost/geometry/strategies/geographic/buffer_end_round.hpp b/include/boost/geometry/strategies/geographic/buffer_end_round.hpp index fbd80c43f..f4b555122 100644 --- a/include/boost/geometry/strategies/geographic/buffer_end_round.hpp +++ b/include/boost/geometry/strategies/geographic/buffer_end_round.hpp @@ -15,6 +15,7 @@ #include #include +#include #include #include #include @@ -34,24 +35,10 @@ template > class geographic_end_round { - static bool const enable_azimuth = true; - static bool const enable_coordinates = true; - - template - using inverse = typename FormulaPolicy::template inverse - < - T, false, enable_azimuth, false, false, false - >; - template - using direct = typename FormulaPolicy::template direct - < - T, enable_coordinates, false, false, false - >; - public : //! \brief Constructs the strategy - //! \param points_per_circle points which would be used for a full circle + //! \param points_per_circle Number of points which would be used for a full circle //! (if points_per_circle is smaller than 4, it is internally set to 4) explicit inline geographic_end_round(std::size_t points_per_circle = 90) : m_points_per_circle((points_per_circle < 4u) ? 4u : points_per_circle) @@ -61,7 +48,7 @@ public : template inline void generate(T lon_rad, T lat_rad, T distance, T azimuth, RangeOut& range_out) const { - using point_t = typename boost::range_value::type; + using helper = geographic_buffer_helper; std::size_t const n = m_points_per_circle / 2; T const angle_diff = geometry::math::pi() / n; T azi = math::wrap_azimuth_in_radian(azimuth + angle_diff); @@ -70,11 +57,7 @@ public : // because left and right are inserted before and after this range. for (std::size_t i = 1; i < n; i++) { - auto const d = direct::apply(lon_rad, lat_rad, distance, azi, m_spheroid); - point_t point; - set_from_radian<0>(point, d.lon2); - set_from_radian<1>(point, d.lat2); - range_out.emplace_back(point); + helper::append_point(lon_rad, lat_rad, distance, azi, m_spheroid, range_out); azi = math::wrap_azimuth_in_radian(azi + angle_diff); } } @@ -93,6 +76,13 @@ public : CalculationType >::type; + using helper = geographic_buffer_helper; + + calc_t const lon_rad = get_as_radian<0>(ultimate_point); + calc_t const lat_rad = get_as_radian<1>(ultimate_point); + + auto const azimuth = helper::azimuth(lon_rad, lat_rad, perp_left_point, m_spheroid); + calc_t const dist_left = distance.apply(penultimate_point, ultimate_point, buffer_side_left); calc_t const dist_right = distance.apply(penultimate_point, ultimate_point, buffer_side_right); @@ -100,13 +90,6 @@ public : || (side == buffer_side_right && dist_left < 0 && -dist_left > dist_right) ; - calc_t const lon_rad = get_as_radian<0>(ultimate_point); - calc_t const lat_rad = get_as_radian<1>(ultimate_point); - calc_t const lon1_rad = get_as_radian<0>(perp_left_point); - calc_t const lat1_rad = get_as_radian<1>(perp_left_point); - - auto const azimuth = inverse::apply(lon_rad, lat_rad, lon1_rad, lat1_rad, m_spheroid).azimuth; - if (reversed) { range_out.push_back(perp_right_point); @@ -129,7 +112,7 @@ public : = (side == buffer_side_right ? (dist_right - dist_left) : (dist_left - dist_right)) / two; - auto const shifted = direct::apply(lon_rad, lat_rad, dist_half, azimuth, m_spheroid); + auto const shifted = helper::direct::apply(lon_rad, lat_rad, dist_half, azimuth, m_spheroid); generate(shifted.lon2, shifted.lat2, dist_average, azimuth, range_out); } @@ -155,7 +138,6 @@ private : Spheroid m_spheroid; }; - }} // namespace strategy::buffer }} // namespace boost::geometry diff --git a/include/boost/geometry/strategies/geographic/buffer_helper.hpp b/include/boost/geometry/strategies/geographic/buffer_helper.hpp new file mode 100644 index 000000000..16caa3d02 --- /dev/null +++ b/include/boost/geometry/strategies/geographic/buffer_helper.hpp @@ -0,0 +1,101 @@ +// Boost.Geometry + +// Copyright (c) 2022 Barend Gehrels, Amsterdam, the Netherlands. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_BUFFER_HELPER_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_BUFFER_HELPER_HPP + +#include +#include +#include +#include + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace buffer +{ + +#ifndef DOXYGEN_SHOULD_SKIP_THIS +template +struct geographic_buffer_helper +{ + static bool const enable_azimuth = true; + static bool const enable_coordinates = true; + + using inverse = typename FormulaPolicy::template inverse + < + CalculationType, false, enable_azimuth, false, false, false + >; + + using direct = typename FormulaPolicy::template direct + < + CalculationType, enable_coordinates, false, false, false + >; + + // Calculates the azimuth using the inverse formula, where the first point + // is specified by lon/lat (for pragmatic reasons) and the second point as a point. + template + static inline CalculationType azimuth(T const& lon_rad, T const& lat_rad, + Point const& p, Spheroid const& spheroid) + { + return inverse::apply(lon_rad, lat_rad, get_as_radian<0>(p), get_as_radian<1>(p), spheroid).azimuth; + } + + // Using specified points, distance and azimuth it calculates a new point + // and appends it to the range + template + static inline void append_point(T const& lon_rad, T const& lat_rad, + T const& distance, T const& angle, + Spheroid const& spheroid, RangeOut& range_out) + { + using point_t = typename boost::range_value::type; + point_t point; + auto const d = direct::apply(lon_rad, lat_rad, distance, angle, spheroid); + set_from_radian<0>(point, d.lon2); + set_from_radian<1>(point, d.lat2); + range_out.emplace_back(point); + } + + // Calculates the angle diff and azimuth of a point (specified as lon/lat) + // and two points, perpendicular in the buffer context. + template + static inline bool calculate_angles(T const& lon_rad, T const& lat_rad, Point const& perp1, + Point const& perp2, Spheroid const& spheroid, + T& angle_diff, T& first_azimuth) + { + T const inv1 = azimuth(lon_rad, lat_rad, perp1, spheroid); + T const inv2 = azimuth(lon_rad, lat_rad, perp2, spheroid); + + static CalculationType const two_pi = geometry::math::two_pi(); + static CalculationType const pi = geometry::math::pi(); + + // For a sharp corner, perpendicular points are nearly opposite and the + // angle between the two azimuths can be nearly 180, but not more. + angle_diff = inv2 < inv1 ? (two_pi + inv2) - inv1 : inv2 - inv1; + + if (angle_diff < 0 || angle_diff > pi) + { + // Defensive check with asserts + BOOST_GEOMETRY_ASSERT(angle_diff >= 0); + BOOST_GEOMETRY_ASSERT(angle_diff <= pi); + return false; + } + + first_azimuth = inv1; + + return true; + } +}; +#endif // DOXYGEN_SHOULD_SKIP_THIS + +}} // namespace strategy::buffer + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_BUFFER_HELPER_HPP diff --git a/include/boost/geometry/strategies/geographic/buffer_join_miter.hpp b/include/boost/geometry/strategies/geographic/buffer_join_miter.hpp index d15214d28..ae4895d0d 100644 --- a/include/boost/geometry/strategies/geographic/buffer_join_miter.hpp +++ b/include/boost/geometry/strategies/geographic/buffer_join_miter.hpp @@ -15,6 +15,7 @@ #include #include +#include #include #include #include @@ -34,20 +35,6 @@ template > class geographic_join_miter { - static bool const enable_azimuth = true; - static bool const enable_coordinates = true; - - template - using inverse = typename FormulaPolicy::template inverse - < - T, false, enable_azimuth, false, false, false - >; - template - using direct = typename FormulaPolicy::template direct - < - T, enable_coordinates, false, false, false - >; - public : //! \brief Constructs the strategy @@ -71,38 +58,24 @@ public : CalculationType >::type; + using helper = geographic_buffer_helper; + calc_t const lon_rad = get_as_radian<0>(vertex); calc_t const lat_rad = get_as_radian<1>(vertex); - calc_t const lon1_rad = get_as_radian<0>(perp1); - calc_t const lat1_rad = get_as_radian<1>(perp1); - calc_t const lon2_rad = get_as_radian<0>(perp2); - calc_t const lat2_rad = get_as_radian<1>(perp2); - // Calculate angles from vertex to perp1/perp2 - auto const inv1 = inverse::apply(lon_rad, lat_rad, lon1_rad, lat1_rad, m_spheroid); - auto const inv2 = inverse::apply(lon_rad, lat_rad, lon2_rad, lat2_rad, m_spheroid); - - // For a sharp corner, perpendicular points are nearly opposite and the - // angle between the two azimuths can be nearly 180, but not more. - calc_t const two_pi = geometry::math::two_pi(); - bool const wrapped = inv2.azimuth < inv1.azimuth; - calc_t const angle_diff = wrapped - ? ((two_pi + inv2.azimuth) - inv1.azimuth) - : inv2.azimuth - inv1.azimuth; - - if (angle_diff < 0 || angle_diff > geometry::math::pi()) + calc_t first_azimuth; + calc_t angle_diff; + if (! helper::calculate_angles(lon_rad, lat_rad, perp1, perp2, m_spheroid, + angle_diff, first_azimuth)) { - // Defensive check with asserts - BOOST_GEOMETRY_ASSERT(angle_diff >= 0); - BOOST_GEOMETRY_ASSERT(angle_diff <= geometry::math::pi()); return false; } calc_t const half = 0.5; calc_t const half_angle_diff = half * angle_diff; - calc_t const cos_angle = std::cos(half_angle_diff); + calc_t const azi = math::wrap_azimuth_in_radian(first_azimuth + half_angle_diff); - calc_t const max_distance = m_miter_limit * geometry::math::abs(buffer_distance); + calc_t const cos_angle = std::cos(half_angle_diff); if (cos_angle == 0) { @@ -111,16 +84,11 @@ public : } // If it is sharp (angle close to 0), the distance will become too high and will be capped. + calc_t const max_distance = m_miter_limit * geometry::math::abs(buffer_distance); calc_t const distance = (std::min)(max_distance, buffer_distance / cos_angle); - calc_t const azi = math::wrap_azimuth_in_radian(inv1.azimuth + half_angle_diff); - - Point point; - auto const d = direct::apply(lon_rad, lat_rad, distance, azi, m_spheroid); - set_from_radian<0>(point, d.lon2); - set_from_radian<1>(point, d.lat2); range_out.push_back(perp1); - range_out.push_back(point); + helper::append_point(lon_rad, lat_rad, distance, azi, m_spheroid, range_out); range_out.push_back(perp2); return true; } diff --git a/include/boost/geometry/strategies/geographic/buffer_join_round.hpp b/include/boost/geometry/strategies/geographic/buffer_join_round.hpp index 1fcf5e7ac..87ae8db97 100644 --- a/include/boost/geometry/strategies/geographic/buffer_join_round.hpp +++ b/include/boost/geometry/strategies/geographic/buffer_join_round.hpp @@ -15,6 +15,7 @@ #include #include +#include #include #include #include @@ -34,26 +35,13 @@ template > class geographic_join_round { - static bool const enable_azimuth = true; - static bool const enable_coordinates = true; - - template - using inverse = typename FormulaPolicy::template inverse - < - T, false, enable_azimuth, false, false, false - >; - template - using direct = typename FormulaPolicy::template direct - < - T, enable_coordinates, false, false, false - >; - public : //! \brief Constructs the strategy - //! \param points_per_circle points which would be used for a full circle + //! \param points_per_circle Number of points which would be used for a full circle + //! (if points_per_circle is smaller than 4, it is internally set to 4) explicit inline geographic_join_round(std::size_t points_per_circle = 90) - : m_points_per_circle(points_per_circle) + : m_points_per_circle((points_per_circle < 4u) ? 4u : points_per_circle) {} #ifndef DOXYGEN_SHOULD_SKIP_THIS @@ -71,39 +59,26 @@ public : CalculationType >::type; + using helper = geographic_buffer_helper; + calc_t const lon_rad = get_as_radian<0>(vertex); calc_t const lat_rad = get_as_radian<1>(vertex); - calc_t const lon1_rad = get_as_radian<0>(perp1); - calc_t const lat1_rad = get_as_radian<1>(perp1); - calc_t const lon2_rad = get_as_radian<0>(perp2); - calc_t const lat2_rad = get_as_radian<1>(perp2); - // Calculate angles from vertex to perp1/perp2 - auto const inv1 = inverse::apply(lon_rad, lat_rad, lon1_rad, lat1_rad, m_spheroid); - auto const inv2 = inverse::apply(lon_rad, lat_rad, lon2_rad, lat2_rad, m_spheroid); - - // For a sharp corner, perpendicular points are nearly opposite and the - // angle between the two azimuths can be nearly 180, but not more. - calc_t const two_pi = geometry::math::two_pi(); - bool const wrapped = inv2.azimuth < inv1.azimuth; - calc_t const angle_diff = wrapped - ? ((two_pi + inv2.azimuth) - inv1.azimuth) - : inv2.azimuth - inv1.azimuth; - - if (angle_diff < 0 || angle_diff > geometry::math::pi()) + calc_t first_azimuth; + calc_t angle_diff; + if (! helper::calculate_angles(lon_rad, lat_rad, perp1, perp2, m_spheroid, + angle_diff, first_azimuth)) { - // Defensive check with asserts - BOOST_GEOMETRY_ASSERT(angle_diff >= 0); - BOOST_GEOMETRY_ASSERT(angle_diff <= geometry::math::pi()); return false; } + static calc_t const two_pi = geometry::math::two_pi(); calc_t const circle_fraction = angle_diff / two_pi; std::size_t const n = (std::max)(static_cast( std::ceil(m_points_per_circle * circle_fraction)), std::size_t(1)); calc_t const diff = angle_diff / static_cast(n); - calc_t azi = math::wrap_azimuth_in_radian(inv1.azimuth + diff); + calc_t azi = math::wrap_azimuth_in_radian(first_azimuth + diff); range_out.push_back(perp1); @@ -111,12 +86,7 @@ public : // because perp1 and perp2 are inserted before and after this range. for (std::size_t i = 1; i < n; i++) { - auto const d = direct::apply(lon_rad, lat_rad, buffer_distance, azi, m_spheroid); - Point p; - set_from_radian<0>(p, d.lon2); - set_from_radian<1>(p, d.lat2); - range_out.emplace_back(p); - + helper::append_point(lon_rad, lat_rad, buffer_distance, azi, m_spheroid, range_out); azi = math::wrap_azimuth_in_radian(azi + diff); } diff --git a/include/boost/geometry/strategies/geographic/buffer_point_circle.hpp b/include/boost/geometry/strategies/geographic/buffer_point_circle.hpp index 85dbe635b..c2d98ea68 100644 --- a/include/boost/geometry/strategies/geographic/buffer_point_circle.hpp +++ b/include/boost/geometry/strategies/geographic/buffer_point_circle.hpp @@ -17,12 +17,16 @@ #include +#include + #include #include +#include #include #include #include + namespace boost { namespace geometry { @@ -55,77 +59,67 @@ template class geographic_point_circle { public : + //! \brief Constructs the strategy - //! \param count number of points for the created circle (if count - //! is smaller than 3, count is internally set to 3) - explicit geographic_point_circle(std::size_t count = 90) - : m_count((count < 3u) ? 3u : count) + //! \param points_per_circle Number of points for a full circle + //! (if points_per_circle is smaller than 3, it is internally set to 3) + explicit geographic_point_circle(std::size_t points_per_circle = 90) + : m_points_per_circle((points_per_circle < 3u) ? 3u : points_per_circle) {} #ifndef DOXYGEN_SHOULD_SKIP_THIS - //! Fills output_range with a circle around point using distance_strategy + //! Fills range_out with a circle around point using distance_strategy template < typename Point, - typename OutputRange, + typename RangeOut, typename DistanceStrategy > inline void apply(Point const& point, DistanceStrategy const& distance_strategy, - OutputRange& output_range) const + RangeOut& range_out) const { - using output_point_type = typename boost::range_value::type; - - using calculation_type = typename select_calculation_type + using calc_t = typename select_calculation_type < - Point, output_point_type, + Point, + typename boost::range_value::type, CalculationType >::type; - auto const lon_rad = get_as_radian<0>(point); - auto const lat_rad = get_as_radian<1>(point); + using helper = geographic_buffer_helper; - calculation_type const buffer_distance = distance_strategy.apply(point, + calc_t const lon_rad = get_as_radian<0>(point); + calc_t const lat_rad = get_as_radian<1>(point); + + calc_t const buffer_distance = distance_strategy.apply(point, point, strategy::buffer::buffer_side_left); - using direct_t = typename FormulaPolicy::template direct - < - calculation_type, true, false, false, false - >; + calc_t const two_pi = geometry::math::two_pi(); + calc_t const pi = geometry::math::pi(); - calculation_type const two_pi = geometry::math::two_pi(); - calculation_type const pi = geometry::math::pi(); + calc_t const diff = two_pi / calc_t(m_points_per_circle); + calc_t angle = -pi; - calculation_type const diff = two_pi / calculation_type(m_count); - calculation_type angle = -pi; - - for (std::size_t i = 0; i < m_count; i++, angle += diff) + for (std::size_t i = 0; i < m_points_per_circle; i++, angle += diff) { // If angle is zero, shift angle a tiny bit to avoid spikes. - calculation_type const eps = angle == 0 ? 1.0e-10 : 0.0; - auto const dir_rad = direct_t::apply(lon_rad, lat_rad, - buffer_distance, angle + eps, - m_spheroid); - output_point_type p; - set_from_radian<0>(p, dir_rad.lon2); - set_from_radian<1>(p, dir_rad.lat2); - output_range.push_back(p); + calc_t const eps = angle == 0 ? 1.0e-10 : 0.0; + helper::append_point(lon_rad, lat_rad, buffer_distance, angle + eps, m_spheroid, range_out); } { // Close the range - auto const p = output_range.front(); - output_range.push_back(p); + auto const p = range_out.front(); + range_out.push_back(p); } } #endif // DOXYGEN_SHOULD_SKIP_THIS private : - std::size_t m_count; + std::size_t m_points_per_circle; Spheroid m_spheroid; }; - }} // namespace strategy::buffer }} // namespace boost::geometry diff --git a/include/boost/geometry/strategies/geographic/buffer_side_straight.hpp b/include/boost/geometry/strategies/geographic/buffer_side_straight.hpp index 37625dc5f..93588e2af 100644 --- a/include/boost/geometry/strategies/geographic/buffer_side_straight.hpp +++ b/include/boost/geometry/strategies/geographic/buffer_side_straight.hpp @@ -13,12 +13,14 @@ #include +#include + #include #include +#include #include #include #include -#include namespace boost { namespace geometry @@ -52,21 +54,23 @@ public : template < typename Point, - typename OutputRange, + typename RangeOut, typename DistanceStrategy > inline result_code apply(Point const& input_p1, Point const& input_p2, buffer_side_selector side, DistanceStrategy const& distance_strategy, - OutputRange& output_range) const + RangeOut& range_out) const { using calc_t = typename select_calculation_type < Point, - typename boost::range_value::type, + typename boost::range_value::type, CalculationType >::type; + using helper = geographic_buffer_helper; + calc_t const lon1_rad = get_as_radian<0>(input_p1); calc_t const lat1_rad = get_as_radian<1>(input_p1); calc_t const lon2_rad = get_as_radian<0>(input_p2); @@ -79,44 +83,16 @@ public : return result_no_output; } - // Define the types for the Formulas to calculate a point - // at a certain distance (using with coordinates) - // and to calculate the angle between two specified points - // (using with azimuth) - // See also boost/geometry/strategies/geographic/parameters.hpp - constexpr bool enable_azimuth = true; - constexpr bool enable_coordinates = true; - - using direct_t = typename FormulaPolicy::template direct - < - calc_t, enable_coordinates, false, false, false - >; - - using inverse_t = typename FormulaPolicy::template inverse - < - calc_t, false, enable_azimuth, false, false, false - >; - // Measure the angle from p1 to p2 with the Inverse transformation, // and subtract pi/2 to make it perpendicular. - auto const inv = inverse_t::apply(lon1_rad, lat1_rad, - lon2_rad, lat2_rad, m_spheroid); - auto const angle = math::wrap_azimuth_in_radian(inv.azimuth - - geometry::math::half_pi()); + auto const inv = helper::azimuth(lon1_rad, lat1_rad, input_p2, m_spheroid); + auto const angle = math::wrap_azimuth_in_radian(inv - geometry::math::half_pi()); // Calculate the distance and generate two points at that distance - // with the Direct transformation auto const distance = distance_strategy.apply(input_p1, input_p2, side); - auto const d1 = direct_t::apply(lon1_rad, lat1_rad, distance, angle, m_spheroid); - auto const d2 = direct_t::apply(lon2_rad, lat2_rad, distance, angle, m_spheroid); - - output_range.resize(2); - - set_from_radian<0>(output_range.front(), d1.lon2); - set_from_radian<1>(output_range.front(), d1.lat2); - set_from_radian<0>(output_range.back(), d2.lon2); - set_from_radian<1>(output_range.back(), d2.lat2); + helper::append_point(lon1_rad, lat1_rad, distance, angle, m_spheroid, range_out); + helper::append_point(lon2_rad, lat2_rad, distance, angle, m_spheroid, range_out); return result_normal; } @@ -126,7 +102,6 @@ private : Spheroid m_spheroid; }; - }} // namespace strategy::buffer }} // namespace boost::geometry