diff options
Diffstat (limited to 'boost/geometry/formulas/area_formulas.hpp')
| -rw-r--r-- | boost/geometry/formulas/area_formulas.hpp | 170 |
1 files changed, 98 insertions, 72 deletions
diff --git a/boost/geometry/formulas/area_formulas.hpp b/boost/geometry/formulas/area_formulas.hpp index d459528c27..b8da8a6727 100644 --- a/boost/geometry/formulas/area_formulas.hpp +++ b/boost/geometry/formulas/area_formulas.hpp @@ -1,7 +1,8 @@ // Boost.Geometry -// Copyright (c) 2015-2021 Oracle and/or its affiliates. +// Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland. +// Copyright (c) 2015-2022 Oracle and/or its affiliates. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -342,6 +343,15 @@ public: } } + static inline CT trapezoidal_formula(CT lat1r, CT lat2r, CT lon21r) + { + CT const c1 = CT(1); + CT const c2 = CT(2); + CT const tan_lat1 = tan(lat1r / c2); + CT const tan_lat2 = tan(lat2r / c2); + + return c2 * atan(((tan_lat1 + tan_lat2) / (c1 + tan_lat1 * tan_lat2))* tan(lon21r / c2)); + } /* Compute the spherical excess of a geodesic (or shperical) segment @@ -354,43 +364,44 @@ public: static inline CT spherical(PointOfSegment const& p1, PointOfSegment const& p2) { + CT const pi = math::pi<CT>(); + CT excess; - if (LongSegment) // not for segments parallel to equator + CT const lon1r = get_as_radian<0>(p1); + CT const lat1r = get_as_radian<1>(p1); + CT const lon2r = get_as_radian<0>(p2); + CT const lat2r = get_as_radian<1>(p2); + + CT lon12r = lon2r - lon1r; + math::normalize_longitude<radian, CT>(lon12r); + + if (lon12r == pi || lon12r == -pi) { - CT cbet1 = cos(geometry::get_as_radian<1>(p1)); - CT sbet1 = sin(geometry::get_as_radian<1>(p1)); - CT cbet2 = cos(geometry::get_as_radian<1>(p2)); - CT sbet2 = sin(geometry::get_as_radian<1>(p2)); + return pi; + } - CT omg12 = geometry::get_as_radian<0>(p1) - - geometry::get_as_radian<0>(p2); - CT comg12 = cos(omg12); - CT somg12 = sin(omg12); + if (BOOST_GEOMETRY_CONDITION(LongSegment) && lat1r != lat2r) // not for segments parallel to equator + { + CT const cbet1 = cos(lat1r); + CT const sbet1 = sin(lat1r); + CT const cbet2 = cos(lat2r); + CT const sbet2 = sin(lat2r); - CT alp1 = atan2(cbet1 * sbet2 - - sbet1 * cbet2 * comg12, - cbet2 * somg12); + CT const omg12 = lon2r - lon1r; + CT const comg12 = cos(omg12); + CT const somg12 = sin(omg12); - CT alp2 = atan2(cbet1 * sbet2 * comg12 - - sbet1 * cbet2, - cbet1 * somg12); + CT const cbet1_sbet2 = cbet1 * sbet2; + CT const sbet1_cbet2 = sbet1 * cbet2; + CT const alp1 = atan2(cbet1_sbet2 - sbet1_cbet2 * comg12, cbet2 * somg12); + CT const alp2 = atan2(cbet1_sbet2 * comg12 - sbet1_cbet2, cbet1 * somg12); excess = alp2 - alp1; } else { - // Trapezoidal formula - - CT tan_lat1 = - tan(geometry::get_as_radian<1>(p1) / 2.0); - CT tan_lat2 = - tan(geometry::get_as_radian<1>(p2) / 2.0); - - excess = CT(2.0) - * atan(((tan_lat1 + tan_lat2) / (CT(1) + tan_lat1 * tan_lat2)) - * tan((geometry::get_as_radian<0>(p2) - - geometry::get_as_radian<0>(p1)) / 2)); + excess = trapezoidal_formula(lat1r, lat2r, lon12r); } return excess; @@ -437,8 +448,13 @@ public: // Constants + CT const c0 = CT(0); + CT const c1 = CT(1); + CT const c2 = CT(2); + CT const pi = math::pi<CT>(); + CT const half_pi = pi / c2; CT const ep = spheroid_const.m_ep; - CT const one_minus_f = CT(1) - spheroid_const.m_f; + CT const one_minus_f = c1 - spheroid_const.m_f; // Basic trigonometric computations // the compiler could optimize here using sincos function @@ -472,43 +488,47 @@ public: CT excess; - auto const half_pi = math::pi<CT>() / 2; - bool meridian = lon2r - lon1r == CT(0) - || lat1r == half_pi || lat1r == -half_pi - || lat2r == half_pi || lat2r == -half_pi; + CT lon12r = lon2r - lon1r; + math::normalize_longitude<radian, CT>(lon12r); - if (!meridian && (i_res.distance) - < mean_radius<CT>(spheroid_const.m_spheroid) / CT(638)) // short segment + // Comparing with "==" works with all test cases here, but could potential create numerical issues + if (lon12r == pi || lon12r == -pi) { - CT tan_lat1 = tan(lat1r / 2.0); - CT tan_lat2 = tan(lat2r / 2.0); - - excess = CT(2.0) - * atan(((tan_lat1 + tan_lat2) / (CT(1) + tan_lat1 * tan_lat2)) - * tan((lon2r - lon1r) / 2)); + result.spherical_term = pi; } else { - /* in some cases this formula gives more accurate results - * - * CT sin_omg12 = cos_omg1 * sin_omg2 - sin_omg1 * cos_omg2; - normalize(sin_omg12, cos_omg12); - - CT cos_omg12p1 = CT(1) + cos_omg12; - CT cos_bet1p1 = CT(1) + cos_bet1; - CT cos_bet2p1 = CT(1) + cos_bet2; - excess = CT(2) * atan2(sin_omg12 * (sin_bet1 * cos_bet2p1 + sin_bet2 * cos_bet1p1), - cos_omg12p1 * (sin_bet1 * sin_bet2 + cos_bet1p1 * cos_bet2p1)); - */ - - excess = alp2 - alp1; + bool const meridian = lon12r == c0 + || lat1r == half_pi || lat1r == -half_pi + || lat2r == half_pi || lat2r == -half_pi; + + if (!meridian && (i_res.distance) + < mean_radius<CT>(spheroid_const.m_spheroid) / CT(638)) // short segment + { + excess = trapezoidal_formula(lat1r, lat2r, lon12r); + } + else + { + /* in some cases this formula gives more accurate results + CT sin_omg12 = cos_omg1 * sin_omg2 - sin_omg1 * cos_omg2; + normalize(sin_omg12, cos_omg12); + + CT cos_omg12p1 = CT(1) + cos_omg12; + CT cos_bet1p1 = CT(1) + cos_bet1; + CT cos_bet2p1 = CT(1) + cos_bet2; + excess = CT(2) * atan2(sin_omg12 * (sin_bet1 * cos_bet2p1 + sin_bet2 * cos_bet1p1), + cos_omg12p1 * (sin_bet1 * sin_bet2 + cos_bet1p1 * cos_bet2p1)); + */ + + excess = alp2 - alp1; + } + + result.spherical_term = excess; } - result.spherical_term = excess; - // Ellipsoidal term computation (uses integral approximation) - CT const cos_alp0 = math::sqrt(CT(1) - math::sqr(sin_alp0)); + CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0)); //CT const cos_alp0 = hypot(cos_alp1, sin_alp1 * sin_bet1); CT cos_sig1 = cos_alp1 * cos_bet1; CT cos_sig2 = cos_alp2 * cos_bet2; @@ -523,8 +543,8 @@ public: if (ExpandEpsN) // expand by eps and n { CT const k2 = math::sqr(ep * cos_alp0); - CT const sqrt_k2_plus_one = math::sqrt(CT(1) + k2); - CT const eps = (sqrt_k2_plus_one - CT(1)) / (sqrt_k2_plus_one + CT(1)); + CT const sqrt_k2_plus_one = math::sqrt(c1 + k2); + CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1); // Generate and evaluate the polynomials on eps (i.e. var2 = eps) // to get the final series coefficients @@ -563,20 +583,26 @@ public: static inline bool crosses_prime_meridian(PointOfSegment const& p1, PointOfSegment const& p2) { - CT const pi - = geometry::math::pi<CT>(); - CT const two_pi - = geometry::math::two_pi<CT>(); - - CT p1_lon = get_as_radian<0>(p1) - - ( floor( get_as_radian<0>(p1) / two_pi ) - * two_pi ); - CT p2_lon = get_as_radian<0>(p2) - - ( floor( get_as_radian<0>(p2) / two_pi ) - * two_pi ); - - CT max_lon = (std::max)(p1_lon, p2_lon); - CT min_lon = (std::min)(p1_lon, p2_lon); + CT const pi = geometry::math::pi<CT>(); + CT const two_pi = geometry::math::two_pi<CT>(); + + CT const lon1r = get_as_radian<0>(p1); + CT const lon2r = get_as_radian<0>(p2); + + CT lon12 = lon2r - lon1r; + math::normalize_longitude<radian, CT>(lon12); + + // Comparing with "==" works with all test cases here, but could potential create numerical issues + if (lon12 == pi || lon12 == -pi) + { + return true; + } + + CT const p1_lon = lon1r - ( std::floor( lon1r / two_pi ) * two_pi ); + CT const p2_lon = lon2r - ( std::floor( lon2r / two_pi ) * two_pi ); + + CT const max_lon = (std::max)(p1_lon, p2_lon); + CT const min_lon = (std::min)(p1_lon, p2_lon); return max_lon > pi && min_lon < pi && max_lon - min_lon > pi; } |