[buffer] harmonize geographic strategy code and share code

This commit is contained in:
Barend Gehrels 2022-10-19 11:43:48 +02:00
parent 4a255ab4b3
commit 1df1e7021e
6 changed files with 179 additions and 189 deletions

View File

@ -15,6 +15,7 @@
#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/strategies/buffer.hpp> #include <boost/geometry/strategies/buffer.hpp>
#include <boost/geometry/strategies/geographic/buffer_helper.hpp>
#include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp>
#include <boost/geometry/util/math.hpp> #include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_calculation_type.hpp> #include <boost/geometry/util/select_calculation_type.hpp>
@ -34,24 +35,10 @@ template
> >
class geographic_end_round class geographic_end_round
{ {
static bool const enable_azimuth = true;
static bool const enable_coordinates = true;
template <typename T>
using inverse = typename FormulaPolicy::template inverse
<
T, false, enable_azimuth, false, false, false
>;
template <typename T>
using direct = typename FormulaPolicy::template direct
<
T, enable_coordinates, false, false, false
>;
public : public :
//! \brief Constructs the strategy //! \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) //! (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) explicit inline geographic_end_round(std::size_t points_per_circle = 90)
: m_points_per_circle((points_per_circle < 4u) ? 4u : points_per_circle) : m_points_per_circle((points_per_circle < 4u) ? 4u : points_per_circle)
@ -61,7 +48,7 @@ public :
template <typename T, typename RangeOut> template <typename T, typename RangeOut>
inline void generate(T lon_rad, T lat_rad, T distance, T azimuth, RangeOut& range_out) const inline void generate(T lon_rad, T lat_rad, T distance, T azimuth, RangeOut& range_out) const
{ {
using point_t = typename boost::range_value<RangeOut const>::type; using helper = geographic_buffer_helper<FormulaPolicy, T>;
std::size_t const n = m_points_per_circle / 2; std::size_t const n = m_points_per_circle / 2;
T const angle_diff = geometry::math::pi<T>() / n; T const angle_diff = geometry::math::pi<T>() / n;
T azi = math::wrap_azimuth_in_radian(azimuth + angle_diff); 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. // because left and right are inserted before and after this range.
for (std::size_t i = 1; i < n; i++) for (std::size_t i = 1; i < n; i++)
{ {
auto const d = direct<T>::apply(lon_rad, lat_rad, distance, azi, m_spheroid); helper::append_point(lon_rad, lat_rad, distance, azi, m_spheroid, range_out);
point_t point;
set_from_radian<0>(point, d.lon2);
set_from_radian<1>(point, d.lat2);
range_out.emplace_back(point);
azi = math::wrap_azimuth_in_radian(azi + angle_diff); azi = math::wrap_azimuth_in_radian(azi + angle_diff);
} }
} }
@ -93,6 +76,13 @@ public :
CalculationType CalculationType
>::type; >::type;
using helper = geographic_buffer_helper<FormulaPolicy, calc_t>;
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_left = distance.apply(penultimate_point, ultimate_point, buffer_side_left);
calc_t const dist_right = distance.apply(penultimate_point, ultimate_point, buffer_side_right); 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) || (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<calc_t>::apply(lon_rad, lat_rad, lon1_rad, lat1_rad, m_spheroid).azimuth;
if (reversed) if (reversed)
{ {
range_out.push_back(perp_right_point); range_out.push_back(perp_right_point);
@ -129,7 +112,7 @@ public :
= (side == buffer_side_right = (side == buffer_side_right
? (dist_right - dist_left) ? (dist_right - dist_left)
: (dist_left - dist_right)) / two; : (dist_left - dist_right)) / two;
auto const shifted = direct<calc_t>::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); generate(shifted.lon2, shifted.lat2, dist_average, azimuth, range_out);
} }
@ -155,7 +138,6 @@ private :
Spheroid m_spheroid; Spheroid m_spheroid;
}; };
}} // namespace strategy::buffer }} // namespace strategy::buffer
}} // namespace boost::geometry }} // namespace boost::geometry

View File

@ -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 <boost/geometry/core/assert.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/strategies/geographic/parameters.hpp>
#include <boost/geometry/util/math.hpp>
namespace boost { namespace geometry
{
namespace strategy { namespace buffer
{
#ifndef DOXYGEN_SHOULD_SKIP_THIS
template <typename FormulaPolicy, typename CalculationType>
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 <typename T, typename Point, typename Spheroid>
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 <typename T, typename Spheroid, typename RangeOut>
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<RangeOut>::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 <typename T, typename Point, typename Spheroid>
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<CalculationType>();
static CalculationType const pi = geometry::math::pi<CalculationType>();
// 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

View File

@ -15,6 +15,7 @@
#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/strategies/buffer.hpp> #include <boost/geometry/strategies/buffer.hpp>
#include <boost/geometry/strategies/geographic/buffer_helper.hpp>
#include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp>
#include <boost/geometry/util/math.hpp> #include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_calculation_type.hpp> #include <boost/geometry/util/select_calculation_type.hpp>
@ -34,20 +35,6 @@ template
> >
class geographic_join_miter class geographic_join_miter
{ {
static bool const enable_azimuth = true;
static bool const enable_coordinates = true;
template <typename T>
using inverse = typename FormulaPolicy::template inverse
<
T, false, enable_azimuth, false, false, false
>;
template <typename T>
using direct = typename FormulaPolicy::template direct
<
T, enable_coordinates, false, false, false
>;
public : public :
//! \brief Constructs the strategy //! \brief Constructs the strategy
@ -71,38 +58,24 @@ public :
CalculationType CalculationType
>::type; >::type;
using helper = geographic_buffer_helper<FormulaPolicy, calc_t>;
calc_t const lon_rad = get_as_radian<0>(vertex); calc_t const lon_rad = get_as_radian<0>(vertex);
calc_t const lat_rad = get_as_radian<1>(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 calc_t first_azimuth;
auto const inv1 = inverse<calc_t>::apply(lon_rad, lat_rad, lon1_rad, lat1_rad, m_spheroid); calc_t angle_diff;
auto const inv2 = inverse<calc_t>::apply(lon_rad, lat_rad, lon2_rad, lat2_rad, m_spheroid); if (! helper::calculate_angles(lon_rad, lat_rad, perp1, perp2, m_spheroid,
angle_diff, first_azimuth))
// 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<calc_t>();
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>())
{ {
// Defensive check with asserts
BOOST_GEOMETRY_ASSERT(angle_diff >= 0);
BOOST_GEOMETRY_ASSERT(angle_diff <= geometry::math::pi<calc_t>());
return false; return false;
} }
calc_t const half = 0.5; calc_t const half = 0.5;
calc_t const half_angle_diff = half * angle_diff; 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) 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. // 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 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<calc_t>::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(perp1);
range_out.push_back(point); helper::append_point(lon_rad, lat_rad, distance, azi, m_spheroid, range_out);
range_out.push_back(perp2); range_out.push_back(perp2);
return true; return true;
} }

View File

@ -15,6 +15,7 @@
#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/strategies/buffer.hpp> #include <boost/geometry/strategies/buffer.hpp>
#include <boost/geometry/strategies/geographic/buffer_helper.hpp>
#include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp>
#include <boost/geometry/util/math.hpp> #include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_calculation_type.hpp> #include <boost/geometry/util/select_calculation_type.hpp>
@ -34,26 +35,13 @@ template
> >
class geographic_join_round class geographic_join_round
{ {
static bool const enable_azimuth = true;
static bool const enable_coordinates = true;
template <typename T>
using inverse = typename FormulaPolicy::template inverse
<
T, false, enable_azimuth, false, false, false
>;
template <typename T>
using direct = typename FormulaPolicy::template direct
<
T, enable_coordinates, false, false, false
>;
public : public :
//! \brief Constructs the strategy //! \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) 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 #ifndef DOXYGEN_SHOULD_SKIP_THIS
@ -71,39 +59,26 @@ public :
CalculationType CalculationType
>::type; >::type;
using helper = geographic_buffer_helper<FormulaPolicy, calc_t>;
calc_t const lon_rad = get_as_radian<0>(vertex); calc_t const lon_rad = get_as_radian<0>(vertex);
calc_t const lat_rad = get_as_radian<1>(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 calc_t first_azimuth;
auto const inv1 = inverse<calc_t>::apply(lon_rad, lat_rad, lon1_rad, lat1_rad, m_spheroid); calc_t angle_diff;
auto const inv2 = inverse<calc_t>::apply(lon_rad, lat_rad, lon2_rad, lat2_rad, m_spheroid); if (! helper::calculate_angles(lon_rad, lat_rad, perp1, perp2, m_spheroid,
angle_diff, first_azimuth))
// 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<calc_t>();
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>())
{ {
// Defensive check with asserts
BOOST_GEOMETRY_ASSERT(angle_diff >= 0);
BOOST_GEOMETRY_ASSERT(angle_diff <= geometry::math::pi<calc_t>());
return false; return false;
} }
static calc_t const two_pi = geometry::math::two_pi<calc_t>();
calc_t const circle_fraction = angle_diff / two_pi; calc_t const circle_fraction = angle_diff / two_pi;
std::size_t const n = (std::max)(static_cast<std::size_t>( std::size_t const n = (std::max)(static_cast<std::size_t>(
std::ceil(m_points_per_circle * circle_fraction)), std::size_t(1)); std::ceil(m_points_per_circle * circle_fraction)), std::size_t(1));
calc_t const diff = angle_diff / static_cast<calc_t>(n); calc_t const diff = angle_diff / static_cast<calc_t>(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); range_out.push_back(perp1);
@ -111,12 +86,7 @@ public :
// because perp1 and perp2 are inserted before and after this range. // because perp1 and perp2 are inserted before and after this range.
for (std::size_t i = 1; i < n; i++) for (std::size_t i = 1; i < n; i++)
{ {
auto const d = direct<calc_t>::apply(lon_rad, lat_rad, buffer_distance, azi, m_spheroid); helper::append_point(lon_rad, lat_rad, buffer_distance, azi, m_spheroid, range_out);
Point p;
set_from_radian<0>(p, d.lon2);
set_from_radian<1>(p, d.lat2);
range_out.emplace_back(p);
azi = math::wrap_azimuth_in_radian(azi + diff); azi = math::wrap_azimuth_in_radian(azi + diff);
} }

View File

@ -17,12 +17,16 @@
#include <boost/range/value_type.hpp> #include <boost/range/value_type.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/strategies/buffer.hpp> #include <boost/geometry/strategies/buffer.hpp>
#include <boost/geometry/strategies/geographic/buffer_helper.hpp>
#include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp>
#include <boost/geometry/util/math.hpp> #include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_calculation_type.hpp> #include <boost/geometry/util/select_calculation_type.hpp>
namespace boost { namespace geometry namespace boost { namespace geometry
{ {
@ -55,77 +59,67 @@ template
class geographic_point_circle class geographic_point_circle
{ {
public : public :
//! \brief Constructs the strategy //! \brief Constructs the strategy
//! \param count number of points for the created circle (if count //! \param points_per_circle Number of points for a full circle
//! is smaller than 3, count is internally set to 3) //! (if points_per_circle is smaller than 3, it is internally set to 3)
explicit geographic_point_circle(std::size_t count = 90) explicit geographic_point_circle(std::size_t points_per_circle = 90)
: m_count((count < 3u) ? 3u : count) : m_points_per_circle((points_per_circle < 3u) ? 3u : points_per_circle)
{} {}
#ifndef DOXYGEN_SHOULD_SKIP_THIS #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 template
< <
typename Point, typename Point,
typename OutputRange, typename RangeOut,
typename DistanceStrategy typename DistanceStrategy
> >
inline void apply(Point const& point, inline void apply(Point const& point,
DistanceStrategy const& distance_strategy, DistanceStrategy const& distance_strategy,
OutputRange& output_range) const RangeOut& range_out) const
{ {
using output_point_type = typename boost::range_value<OutputRange>::type; using calc_t = typename select_calculation_type
using calculation_type = typename select_calculation_type
< <
Point, output_point_type, Point,
typename boost::range_value<RangeOut>::type,
CalculationType CalculationType
>::type; >::type;
auto const lon_rad = get_as_radian<0>(point); using helper = geographic_buffer_helper<FormulaPolicy, calc_t>;
auto const lat_rad = get_as_radian<1>(point);
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); point, strategy::buffer::buffer_side_left);
using direct_t = typename FormulaPolicy::template direct calc_t const two_pi = geometry::math::two_pi<calc_t>();
< calc_t const pi = geometry::math::pi<calc_t>();
calculation_type, true, false, false, false
>;
calculation_type const two_pi = geometry::math::two_pi<calculation_type>(); calc_t const diff = two_pi / calc_t(m_points_per_circle);
calculation_type const pi = geometry::math::pi<calculation_type>(); calc_t angle = -pi;
calculation_type const diff = two_pi / calculation_type(m_count); for (std::size_t i = 0; i < m_points_per_circle; i++, angle += diff)
calculation_type angle = -pi;
for (std::size_t i = 0; i < m_count; i++, angle += diff)
{ {
// If angle is zero, shift angle a tiny bit to avoid spikes. // If angle is zero, shift angle a tiny bit to avoid spikes.
calculation_type const eps = angle == 0 ? 1.0e-10 : 0.0; calc_t const eps = angle == 0 ? 1.0e-10 : 0.0;
auto const dir_rad = direct_t::apply(lon_rad, lat_rad, helper::append_point(lon_rad, lat_rad, buffer_distance, angle + eps, m_spheroid, range_out);
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);
} }
{ {
// Close the range // Close the range
auto const p = output_range.front(); auto const p = range_out.front();
output_range.push_back(p); range_out.push_back(p);
} }
} }
#endif // DOXYGEN_SHOULD_SKIP_THIS #endif // DOXYGEN_SHOULD_SKIP_THIS
private : private :
std::size_t m_count; std::size_t m_points_per_circle;
Spheroid m_spheroid; Spheroid m_spheroid;
}; };
}} // namespace strategy::buffer }} // namespace strategy::buffer
}} // namespace boost::geometry }} // namespace boost::geometry

View File

@ -13,12 +13,14 @@
#include <boost/range/value_type.hpp> #include <boost/range/value_type.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/strategies/buffer.hpp> #include <boost/geometry/strategies/buffer.hpp>
#include <boost/geometry/strategies/geographic/buffer_helper.hpp>
#include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp>
#include <boost/geometry/util/math.hpp> #include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_calculation_type.hpp> #include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/algorithms/azimuth.hpp>
namespace boost { namespace geometry namespace boost { namespace geometry
@ -52,21 +54,23 @@ public :
template template
< <
typename Point, typename Point,
typename OutputRange, typename RangeOut,
typename DistanceStrategy typename DistanceStrategy
> >
inline result_code apply(Point const& input_p1, Point const& input_p2, inline result_code apply(Point const& input_p1, Point const& input_p2,
buffer_side_selector side, buffer_side_selector side,
DistanceStrategy const& distance_strategy, DistanceStrategy const& distance_strategy,
OutputRange& output_range) const RangeOut& range_out) const
{ {
using calc_t = typename select_calculation_type using calc_t = typename select_calculation_type
< <
Point, Point,
typename boost::range_value<OutputRange>::type, typename boost::range_value<RangeOut>::type,
CalculationType CalculationType
>::type; >::type;
using helper = geographic_buffer_helper<FormulaPolicy, calc_t>;
calc_t const lon1_rad = get_as_radian<0>(input_p1); 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 lat1_rad = get_as_radian<1>(input_p1);
calc_t const lon2_rad = get_as_radian<0>(input_p2); calc_t const lon2_rad = get_as_radian<0>(input_p2);
@ -79,44 +83,16 @@ public :
return result_no_output; return result_no_output;
} }
// Define the types for the Formulas to calculate a point
// at a certain distance (using <direct> with coordinates)
// and to calculate the angle between two specified points
// (using <inverse> 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, // Measure the angle from p1 to p2 with the Inverse transformation,
// and subtract pi/2 to make it perpendicular. // and subtract pi/2 to make it perpendicular.
auto const inv = inverse_t::apply(lon1_rad, lat1_rad, auto const inv = helper::azimuth(lon1_rad, lat1_rad, input_p2, m_spheroid);
lon2_rad, lat2_rad, m_spheroid); auto const angle = math::wrap_azimuth_in_radian(inv - geometry::math::half_pi<calc_t>());
auto const angle = math::wrap_azimuth_in_radian(inv.azimuth -
geometry::math::half_pi<calc_t>());
// Calculate the distance and generate two points at that distance // 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 distance = distance_strategy.apply(input_p1, input_p2, side);
auto const d1 = direct_t::apply(lon1_rad, lat1_rad, distance, angle, m_spheroid); helper::append_point(lon1_rad, lat1_rad, distance, angle, m_spheroid, range_out);
auto const d2 = direct_t::apply(lon2_rad, lat2_rad, distance, angle, m_spheroid); helper::append_point(lon2_rad, lat2_rad, distance, angle, m_spheroid, range_out);
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);
return result_normal; return result_normal;
} }
@ -126,7 +102,6 @@ private :
Spheroid m_spheroid; Spheroid m_spheroid;
}; };
}} // namespace strategy::buffer }} // namespace strategy::buffer
}} // namespace boost::geometry }} // namespace boost::geometry