diff --git a/example/compare_length_of_two_segments.cpp b/example/compare_length_of_two_segments.cpp new file mode 100644 index 0000000000..8b2e9316f4 --- /dev/null +++ b/example/compare_length_of_two_segments.cpp @@ -0,0 +1,267 @@ +/************************************************************************* + > File Name: compare_length_of_two_segments.cpp + > Author: Rylynnn + > Mail: jingry0321@gmail.com + > Created Time: Wed 14 Jun 2017 06:38:20 AM CST + ************************************************************************/ +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include + +#define MAX 1000 +#define EPS 1e-9 + +namespace bg = boost::geometry; + +typedef bg::model::point + + > point_type; + +double const earth_f = 1 / 298.257223563; +double const earth_e2 = earth_f * (2 - earth_f); +double const earth_r = 6317.0; + +/*\brief Compare length of two segments using Local, + * flat earth approximation used in the special case + * when P1=P3(the two segments have a common endpoint) + *\see www.edwilliams.org/avform.htm#flat + * */ +int compare_length_of_two_segments_local_earth(point_type p1, point_type p2, + point_type p4) +{ + double distance_result1; + double distance_result2; + double lon1, lat1, lon2, lat2, lon4, lat4; + double dlat, dlon, R1, R2, dis_North, dis_East; + lon1 = bg::get_as_radian<0>(p1); + lat1 = bg::get_as_radian<1>(p1); + lon2 = bg::get_as_radian<0>(p2); + lat2 = bg::get_as_radian<1>(p2); + lon4 = bg::get_as_radian<0>(p4); + lat4 = bg::get_as_radian<1>(p4); + + dlat = lat2 - lat1; + dlon = lon2 - lon1; + + R1 = earth_r * (1 - earth_e2) / + bg::math::sqrt((1 - earth_e2 * bg::math::sqr(sin(lat1))) + * (1 - earth_e2 * bg::math::sqr(sin(lat1))) + * (1 - earth_e2 * bg::math::sqr(sin(lat1)))); + R2 = earth_r / bg::math::sqrt(1 - earth_e2 * bg::math::sqr(sin(lat1))); + + dis_North = R1 * dlat; + dis_East = R2 * cos(lat1) * dlon; + + distance_result1 = bg::math::sqrt(bg::math::sqr(dis_North) + bg::math::sqr(dis_East)); + + dlat = lat2 - lat1; + dlon = lon2 - lon1; + + R1 = earth_r * (1 - earth_e2) / + bg::math::sqrt((1 - earth_e2 * bg::math::sqr(sin(lat1))) + * (1 - earth_e2 * bg::math::sqr(sin(lat1))) + * (1 - earth_e2 * bg::math::sqr(sin(lat1)))); + R2 = earth_r / bg::math::sqrt(1 - earth_e2 * bg::math::sqr(sin(lat1))); + + dis_North = R1 * dlat; + dis_East = R2 * cos(lat1) * dlon; + + distance_result2 = bg::math::sqrt(bg::math::sqr(dis_North) + bg::math::sqr(dis_East)); + + double result = distance_result1 - distance_result2; + + if (result < -EPS) + { + return 1; + } + else if (result > EPS) + { + return 2; + } + else if (fabs(result) < EPS) + { + return 3; + } + +} +/*\brief Compare length of two segments using haversine + *\see http://en.wikipedia.org/wiki/Haversine_formula + * */ +int compare_length_of_two_segments_haversine(point_type p1, point_type p2, + point_type p3, point_type p4) +{ + double distance_result1; + double distance_result2; + double lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4; + lon1 = bg::get_as_radian<0>(p1); + lat1 = bg::get_as_radian<1>(p1); + lon2 = bg::get_as_radian<0>(p2); + lat2 = bg::get_as_radian<1>(p2); + lon3 = bg::get_as_radian<0>(p3); + lat3 = bg::get_as_radian<1>(p3); + lon4 = bg::get_as_radian<0>(p4); + lat4 = bg::get_as_radian<1>(p4); + distance_result1 = bg::math::hav(lat2 - lat1) + + cos(lat1) * cos(lat2) * bg::math::hav(lon2 - lon1); + distance_result2 = bg::math::hav(lat4 - lat3) + + cos(lat3) * cos(lat4) * bg::math::hav(lon4 - lon3); + + double result = distance_result1 - distance_result2; + + if (result < -EPS) + { + return 1; + } + else if (result > EPS) + { + return 2; + } + else if (fabs(result) < EPS) + { + return 3; + } +} +/*\brief Compare length of two segments using cartesian distance + * */ +int compare_length_of_two_segments_cartesian_distance(point_type p1, point_type p2, + point_type p3, point_type p4) +{ + double distance_result1; + double distance_result2; + + distance_result1 = bg::distance + (p1, p2, bg::strategy::distance::pythagoras()); + distance_result2 = bg::distance + (p3, p4, bg::strategy::distance::pythagoras()); + + double result = distance_result1 - distance_result2; + + if (result < -EPS) + { + return 1; + } + else if (result > EPS) + { + return 2; + } + else if (fabs(result) < EPS) + { + return 3; + } +} +int main() +{ + double lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4; + freopen("testin_cd.in", "r", stdin); + freopen("testout_cd.out", "w", stdout); + while(std::cin >> lon1 >> lat1 >> lon2 >> lat2 >> lon3 >> lat3 >> lon4 >> lat4){ + + point_type P1(lon1, lat1); + point_type P2(lon2, lat2); + point_type P3(lon3, lat3); + point_type P4(lon4, lon4); + + bg::srs::spheroid earth; + + std::clock_t distance_start = std::clock(); + double distance_result1; + double distance_result2; + int distance_result; + + for (int i=0; i >(earth)); + distance_result2 = bg::distance + (P3, P4, bg::strategy::distance::vincenty >(earth)); + } + if (distance_result1 - distance_result2 < -EPS) + { + distance_result = 1; + } + else if (distance_result1 - distance_result2 > EPS) + { + distance_result = 2; + } + else if (fabs(distance_result1 - distance_result2) < EPS) + { + distance_result = 3; + } + + std::clock_t distance_stop = std::clock(); + double secs_distance = double(distance_stop - distance_start) / (double)CLOCKS_PER_SEC; + + int comparable_geographic_distance_result; + double secs_comparable_geographic_distance; + + if (lon1 != lon3 || lat1 != lat3) + { + std::clock_t comparable_geographic_distance_start = std::clock(); + + for (int i=0; i File Name: tc_equator.cpp + > Author: Rylynnn + > Mail: jingry0321@gmail.com + > Created Time: 2017年07月12日 星期三 19时29分46秒 + ************************************************************************/ +#include +#include +#include +#include +#include +using namespace std; +int main() +{ + freopen("tc_begin_equator.in", "w", stdout); + double a, b, c, d, e, f, g, h; + for (int i=0; i<90; i++){ + a = 0; + b = 0; + c = 180; + d = i; + for (int j=0; j<90; j++){ + e = 10; + f = 0; + g = 170; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + for (int j=90; j>0; j--){ + e = 10; + f = 0; + g = 10; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + } + for (int i=90; i>0; i--){ + a = 0; + b = 0; + c = 0; + d = i; + for (int j=0; j<90; j++){ + e = 10; + f = 0; + g = 170; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + for (int j=90; j>0; j--){ + e = 10; + f = 0; + g = 10; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + } + return 0; +} diff --git a/example/test_approximation/test_begin_equator/tc_equator_flat.cpp b/example/test_approximation/test_begin_equator/tc_equator_flat.cpp new file mode 100644 index 0000000000..cdbb3774b8 --- /dev/null +++ b/example/test_approximation/test_begin_equator/tc_equator_flat.cpp @@ -0,0 +1,70 @@ +/************************************************************************* + > File Name: tc_equator_flat.cpp + > Author: Rylynnn + > Mail: jingry0321@gmail.com + > Created Time: 2017年07月13日 星期四 00时55分03秒 + ************************************************************************/ +#include +#include +#include +#include +#include +using namespace std; +int main() +{ + freopen("tc_begin_equator_flat.in", "w", stdout); + double a, b, c, d, e, f, g, h; + for (int i=0; i<90; i++){ + a = 0; + b = 0; + c = 180; + d = i; + for (int j=0; j<90; j++){ + e = 0; + f = 0; + g = 170; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + for (int j=90; j>0; j--){ + e = 0; + f = 0; + g = 10; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + } + for (int i=90; i>0; i--){ + a = 0; + b = 0; + c = 0; + d = i; + for (int j=0; j<90; j++){ + e = 0; + f = 0; + g = 170; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + for (int j=90; j>0; j--){ + e = 0; + f = 0; + g = 10; + h = j; + cout << a << ' ' << b << " " + << c << " " << d << ' ' + << e << ' ' << f << ' ' + << g << ' ' << h << endl; + } + } + return 0; +} diff --git a/example/test_approximation/test_begin_equator/testcase_begin_equator.cpp b/example/test_approximation/test_begin_equator/testcase_begin_equator.cpp new file mode 100644 index 0000000000..28a045f4a2 --- /dev/null +++ b/example/test_approximation/test_begin_equator/testcase_begin_equator.cpp @@ -0,0 +1,121 @@ +/************************************************************************* + > File Name: testcase_begin_equator.cpp + > Author: Rylynnn + > Mail: jingry0321@gmail.com + > Created Time: Wed 14 Jun 2017 06:38:20 AM CST + ************************************************************************/ +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#define MAX 1000 +#define EPS 1e-9 + +namespace bg = boost::geometry; + +typedef bg::model::point + + > point_type; + +double const earth_f = 1 / 298.257223563; +double const earth_e2 = earth_f * (2 - earth_f); +double const earth_r = 6317.0; +int main() +{ + double lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4; + freopen("tc_begin_equator.in", "r", stdin); + freopen("tc_begin_equator.out", "w", stdout); + while(std::cin >> lon1 >> lat1 >> lon2 >> lat2 >> lon3 >> lat3 >> lon4 >> lat4){ + point_type P1(lon1, lat1); + point_type P2(lon2, lat2); + point_type P3(lon3, lat3); + point_type P4(lon4, lon4); + + bg::srs::spheroid earth; + + std::clock_t distance_start = std::clock(); + double distance_result1; + double distance_result2; + int distance_result; + + for (int i=0; i >(earth)); + distance_result2 = bg::distance + (P3, P4, bg::strategy::distance::vincenty >(earth)); + } + if (distance_result1 - distance_result2 < -EPS) + { + distance_result = 1; + } + else if (distance_result1 - distance_result2 > EPS) + { + distance_result = 2; + } + else if (fabs(distance_result1 - distance_result2) < EPS) + { + distance_result = 3; + } + + std::clock_t distance_stop = std::clock(); + double secs_distance = double(distance_stop - distance_start) / (double)CLOCKS_PER_SEC; + + int cg_distance_result ; + double secs_comparable_geographic_distance; + + if (lon1 != lon3 || lat1 != lat3) + { + std::clock_t comparable_geographic_distance_start = std::clock(); + + for (int i=0; i + ::apply(P1, P2, P3, P4, earth); + } + std::clock_t comparable_geographic_distance_stop = std::clock(); + secs_comparable_geographic_distance = double(comparable_geographic_distance_stop - + comparable_geographic_distance_start) + / (double)CLOCKS_PER_SEC; + if(distance_result != cg_distance_result){ + std::cout << "comparable_geographic_distancie(Haversine):" + << cg_distance_result << std::endl; + std::cout << "distance: " << distance_result << std::endl; + std::cout << "P1: " << lon1 << ", " << lat1 + << " P2: " << lon2 << ", " << lat2 << std::endl; + std::cout << "P3: " << lon3 << ", " << lat3 + << " P4: " << lon4 << ", " << lat4 << std::endl; + + std::cout << "consistency:" + << (distance_result == cg_distance_result) + << std::endl; + + std::cout << "distance (sec):" << secs_distance << std::endl; + std::cout << "comparable_geographic_distance (sec):" + << secs_comparable_geographic_distance + << std::endl << std::endl; + } + } + } + return 0; +} diff --git a/example/test_approximation/test_begin_equator/testcase_begin_equator_flat.cpp b/example/test_approximation/test_begin_equator/testcase_begin_equator_flat.cpp new file mode 100644 index 0000000000..0a8eb32e4e --- /dev/null +++ b/example/test_approximation/test_begin_equator/testcase_begin_equator_flat.cpp @@ -0,0 +1,120 @@ +/************************************************************************* + > File Name: testcase_begin_equator_flat.cpp + > Author: Rylynnn + > Mail: jingry0321@gmail.com + > Created Time: 2017年07月13日 星期四 00时51分31秒 + ************************************************************************/ +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#define MAX 1000 +#define EPS 1e-9 + +namespace bg = boost::geometry; + +typedef bg::model::point + + > point_type; + +double const earth_f = 1 / 298.257223563; +double const earth_e2 = earth_f * (2 - earth_f); +double const earth_r = 6317.0; +int main() +{ + double lon1, lat1, lon2, lat2, lon3, lat3, lon4, lat4; + freopen("tc_begin_equator_flat.in", "r", stdin); + freopen("tc_begin_equator_flat.out", "w", stdout); + while(std::cin >> lon1 >> lat1 >> lon2 >> lat2 >> lon3 >> lat3 >> lon4 >> lat4){ + point_type P1(lon1, lat1); + point_type P2(lon2, lat2); + point_type P3(lon3, lat3); + point_type P4(lon4, lon4); + + bg::srs::spheroid earth; + + std::clock_t distance_start = std::clock(); + double distance_result1; + double distance_result2; + int distance_result; + + for (int i=0; i >(earth)); + distance_result2 = bg::distance + (P3, P4, bg::strategy::distance::vincenty >(earth)); + } + if (distance_result1 - distance_result2 < -EPS) + { + distance_result = 1; + } + else if (distance_result1 - distance_result2 > EPS) + { + distance_result = 2; + } + else if (fabs(distance_result1 - distance_result2) < EPS) + { + distance_result = 3; + } + std::clock_t distance_stop = std::clock(); + double secs_distance = double(distance_stop - distance_start) / (double)CLOCKS_PER_SEC; + + int cg_distance_result ; + double secs_comparable_geographic_distance; + + if (lon1 == lon3 && lat1 == lat3) + { + std::clock_t comparable_geographic_distance_start = std::clock(); + + for (int i=0; i::apply(P1, P2, P4, earth); + } + std::clock_t comparable_geographic_distance_stop = std::clock(); + secs_comparable_geographic_distance = double(comparable_geographic_distance_stop - + comparable_geographic_distance_start) + / (double)CLOCKS_PER_SEC; + if(distance_result == cg_distance_result){ + std::cout << "comparable_geographic_distance(Flat_approximation):" + << cg_distance_result << std::endl; + std::cout << "distance: " << distance_result << std::endl; + std::cout << "P1: " << lon1 << ", " << lat1 + << " P2: " << lon2 << ", " << lat2 << std::endl; + std::cout << "P3: " << lon3 << ", " << lat3 + << " P4: " << lon4 << ", " << lat4 << std::endl; + + std::cout << "consistency:" + << (distance_result == cg_distance_result) + << std::endl; + + std::cout << "distance (sec):" << secs_distance << std::endl; + std::cout << "comparable_geographic_distance (sec):" + << secs_comparable_geographic_distance + << std::endl << std::endl; + + } + } + } + return 0; +} diff --git a/example/test_geographic.cpp b/example/test_geographic.cpp new file mode 100644 index 0000000000..d7409f22ef --- /dev/null +++ b/example/test_geographic.cpp @@ -0,0 +1,84 @@ +/************************************************************************* + > File Name: test_geographic.cpp + > Author: Rylynnn + > Mail: jingry0321@gmail.com + > Created Time: Wed 14 Jun 2017 06:38:20 AM CST + ************************************************************************/ +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#define MAX 1000 + +namespace bg = boost::geometry; + +int main() +{ + + typedef bg::model::point + + > point_type; + + point_type P1(0, 0, 0); + point_type P2(10, 10, 10); + point_type P3(0, 0.1, 9); + point_type P4(10, 10.1, 10.0); + + bg::srs::spheroid earth; + + std::clock_t distance_start = std::clock(); + double distance_result1; + double distance_result2; + + for (int i=0; i >(earth)); + distance_result2 = bg::distance + (P3, P4, bg::strategy::distance::vincenty >(earth)); + } + std::clock_t distance_stop = std::clock(); + double secs_distance = double(distance_stop - distance_start) / (double)CLOCKS_PER_SEC; + + + std::clock_t comparable_geographic_distance_start = std::clock(); + double comparable_geographic_distance_result1; + double comparable_geographic_distance_result2; + + for (int i=0; i >(earth)); + comparable_geographic_distance_result2 = bg::comparable_geographic_distance + (P3, P4, bg::strategy::distance::vincenty >(earth)); + } + std::clock_t comparable_geographic_distance_stop = std::clock(); + double secs_comparable_geographic_distance = double(comparable_geographic_distance_stop - comparable_geographic_distance_start) + / (double)CLOCKS_PER_SEC; + + + std::cout << "distance 1:" << distance_result1 << std::endl; + std::cout << "comparable_geographic_distance 1:" << comparable_geographic_distance_result1 << std::endl; + + std::cout << "distance 2:" << distance_result2 << std::endl; + std::cout << "comparable_geographic_distance 2:" << comparable_geographic_distance_result2 << std::endl; + + std::cout << "Consistency:" << ((distance_result1 < distance_result2) + == (comparable_geographic_distance_result1 < comparable_geographic_distance_result2)) + << std::endl; + + std::cout << "distance (sec):" << secs_distance << std::endl; + std::cout << "comparable_geographic_distance (sec):" << secs_comparable_geographic_distance << std::endl; + + return 0; +} diff --git a/include/boost/geometry/algorithms/comparable_geographic_distance.hpp b/include/boost/geometry/algorithms/comparable_geographic_distance.hpp new file mode 100644 index 0000000000..726e4cbe7e --- /dev/null +++ b/include/boost/geometry/algorithms/comparable_geographic_distance.hpp @@ -0,0 +1,107 @@ +//Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2014 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2014 Mateusz Loskot, London, UK. + +// This file was modified by Oracle on 2014. +// Modifications copyright (c) 2014, Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle + +// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library +// (geolib/GGL), copyright (c) 1995-2010 Geodan, 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_ALGORITHMS_COMPARABLE_GEOGRAPHIC_DISTANCE_HPP +#define BOOST_GEOMETRY_ALGORITHMS_COMPARABLE_GEOGRAPHIC_DISTANCE_HPP + +#include + +#include +#include +#include + +#include + + +namespace boost { namespace geometry +{ + +namespace resolve_strategy +{ + +struct comparable_geographic_distance +{ + template + static inline + typename comparable_distance_result::type + apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Strategy const& strategy) + { + typedef typename strategy::distance::services::comparable_type + < + Strategy + >::type comparable_strategy_type; + + return dispatch::distance + < + Geometry1, Geometry2, comparable_strategy_type + >::apply(geometry1, + geometry2, + strategy::distance::services::get_comparable + < + Strategy + >::apply(strategy)); + } +}; + +} // namespace resolve_strategy + + +namespace resolve_variant +{ + + +template +struct comparable_geographic_distance +{ + static inline + typename comparable_distance_result::type + apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Strategy const& strategy) + { + return resolve_strategy::comparable_geographic_distance::apply(geometry1, + geometry2, + strategy); + } +}; + + +} // namespace resolve_variant + + +template +inline typename comparable_distance_result::type +comparable_geographic_distance(Geometry1 const& geometry1, Geometry2 const& geometry2, + Strategy const& strategy) +{ + concepts::check(); + concepts::check(); + + return resolve_variant::comparable_geographic_distance + < + Geometry1, + Geometry2, + Strategy + >::apply(geometry1, geometry2, strategy); +} + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_COMPARABLE_DISTANCE_HPP diff --git a/include/boost/geometry/formulas/cartesine_distance.hpp b/include/boost/geometry/formulas/cartesine_distance.hpp new file mode 100644 index 0000000000..dc6adff835 --- /dev/null +++ b/include/boost/geometry/formulas/cartesine_distance.hpp @@ -0,0 +1,65 @@ +// Boost.Geometry + +// Copyright (c) 2015-2017 Oracle and/or its affiliates. + +// Contributed and/or modified by Ruoyun Jing, on behalf of Oracle + +// 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_FORMULAS_CARTESINE_DISTANCE_HPP +#define BOOST_GEOMETRY_FORMULAS_CARTESINE_DISTANCE_HPP + +#include + +#define BOOST_GEOMETRY_EPS 1e-9 + +namespace bg = boost::geometry; +namespace boost { namespace geometry { namespace formula +{ + +template +class cartesine_distance +{ + +public: + typedef int result_type; + template + < + typename T, + typename Spheriod + > + static inline result_type apply(T const& p1, + T const& p2, + T const& p3, + T const& p4, + Spheriod const& spheriod) + { + result_type result; + + CT const distance_result1 = bg::distance + (p1, p2, bg::strategy::distance::pythagoras()); + CT const distance_result2 = bg::distance + (p3, p4, bg::strategy::distance::pythagoras()); + + CT const sub = distance_result1 - distance_result2; + + if (sub < -BOOST_GEOMETRY_EPS) + { + result = 1; + } + else if (sub > BOOST_GEOMETRY_EPS) + { + result = 2; + } + else if (fabs(sub) < BOOST_GEOMETRY_EPS) + { + result = 3; + } + } +}; + +}}} // namespace boost::geometry::formula + +#endif //BOOST_GEOMETRY_FORMULAS_CARTESINE_DISTANCE_HPP diff --git a/include/boost/geometry/formulas/compare_length_andoyer.hpp b/include/boost/geometry/formulas/compare_length_andoyer.hpp new file mode 100644 index 0000000000..117b5bb190 --- /dev/null +++ b/include/boost/geometry/formulas/compare_length_andoyer.hpp @@ -0,0 +1,144 @@ +// Boost.Geometry + +// Copyright (c) 2015-2017 Oracle and/or its affiliates. + +// Contributed and/or modified by Ruoyun Jing, on behalf of Oracle + +// 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_FORMULAS_COMPARE_LENGTH_ANDOYER_HPP +#define BOOST_GEOMETRY_FORMULAS_COMPARE_LENGTH_ANDOYER_HPP + +#include +#include +#include + +#include +#include + +#include + +#define BOOST_GEOMETRY_EPS 1e-9 + +namespace bg = boost::geometry; +namespace boost { namespace geometry { namespace formula +{ + +template +< + typename CT, + typename Spheriod +> +class compare_length_andoyer +{ + +public: + typedef int result_type; + template + < + typename T + > + static inline result_type apply(T const& p1, + T const& p2, + T const& p3, + T const& p4, + Spheriod const& spheriod) + { + result_type result; + + CT const lon1 = bg::get_as_radian<0>(p1); + CT const lat1 = bg::get_as_radian<1>(p1); + CT const lon2 = bg::get_as_radian<0>(p2); + CT const lat2 = bg::get_as_radian<1>(p2); + CT const lon3 = bg::get_as_radian<0>(p3); + CT const lat3 = bg::get_as_radian<1>(p3); + CT const lon4 = bg::get_as_radian<0>(p4); + CT const lat4 = bg::get_as_radian<1>(p4); + + CT const distance_result1 = andoyer_distance(lon1, lat1, lon2, lat2, spheriod); + CT const distance_result2 = andoyer_distance(lon3, lat3, lon4, lat4, spheriod); + + CT const sub = distance_result1 - distance_result2; + + if (sub < -BOOST_GEOMETRY_EPS) + { + result = 1; + } + else if (sub > BOOST_GEOMETRY_EPS) + { + result = 2; + } + else if (fabs(sub) < BOOST_GEOMETRY_EPS) + { + result = 3; + } + return result; + } + +private: + static inline CT andoyer_distance(CT const& lo1, + CT const& la1, + CT const& lo2, + CT const& la2, + Spheriod const& spher) + { + CT const c0 = CT(0); + + if ( math::equals(lo1, lo2) && math::equals(la1, la2) ) + { + return c0; + } + + CT const c1 = CT(1); + CT const pi = math::pi(); + CT const f = formula::flattening(spher); + + CT const dlon = lo2 - lo1; + CT const sin_dlon = sin(dlon); + CT const cos_dlon = cos(dlon); + CT const sin_lat1 = sin(la1); + CT const cos_lat1 = cos(la1); + CT const sin_lat2 = sin(la2); + CT const cos_lat2 = cos(la2); + + // H,G,T = infinity if cos_d = 1 or cos_d = -1 + // lat1 == +-90 && lat2 == +-90 + // lat1 == lat2 && lon1 == lon2 + CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; + // on some platforms cos_d may be outside valid range + if (cos_d < -c1) + cos_d = -c1; + else if (cos_d > c1) + cos_d = c1; + + CT const d = acos(cos_d); // [0, pi] + CT const sin_d = sin(d); // [-1, 1] + + CT const K = math::sqr(sin_lat1-sin_lat2); + CT const L = math::sqr(sin_lat1+sin_lat2); + CT const three_sin_d = CT(3) * sin_d; + + CT const one_minus_cos_d = c1 - cos_d; + CT const one_plus_cos_d = c1 + cos_d; + // cos_d = 1 or cos_d = -1 means that the points are antipodal + + CT const H = math::equals(one_minus_cos_d, c0) ? + c0 : + (d + three_sin_d) / one_minus_cos_d; + CT const G = math::equals(one_plus_cos_d, c0) ? + c0 : + (d - three_sin_d) / one_plus_cos_d; + + CT const dd = -(f/CT(4))*(H*K+G*L); + + CT const a = get_radius<0>(spher); + + return a * (d + dd); + } +}; + +}}} // namespace boost::geometry::formula + +#endif // BOOST_GEOMETRY_FORMULAS_COMPARE_LENGTH_ANDOYER_HPP diff --git a/include/boost/geometry/formulas/compare_length_haversine.hpp b/include/boost/geometry/formulas/compare_length_haversine.hpp new file mode 100644 index 0000000000..585ae52fbb --- /dev/null +++ b/include/boost/geometry/formulas/compare_length_haversine.hpp @@ -0,0 +1,82 @@ +// Boost.Geometry + +// Copyright (c) 2015-2017 Oracle and/or its affiliates. + +// Contributed and/or modified by Ruoyun Jing, on behalf of Oracle + +// 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_FORMULAS_COMPARE_DISTANCE_HAVERSINE_HPP +#define BOOST_GEOMETRY_FORMULAS_COMPARE_DISTANCE_HAVERSINE_HPP + +#include +#include +#include + +#include +#include + +#include + +#define BOOST_GEOMETRY_EPS 1e-9 + +namespace bg = boost::geometry; +namespace boost { namespace geometry { namespace formula +{ + +template +class compare_length_haversine +{ + +public: + typedef int result_type; + template + < + typename T, + typename Spheriod + > + static inline result_type apply(T const& p1, + T const& p2, + T const& p3, + T const& p4, + Spheriod const& spheriod) + { + result_type result; + + CT const lon1 = bg::get_as_radian<0>(p1); + CT const lat1 = bg::get_as_radian<1>(p1); + CT const lon2 = bg::get_as_radian<0>(p2); + CT const lat2 = bg::get_as_radian<1>(p2); + CT const lon3 = bg::get_as_radian<0>(p3); + CT const lat3 = bg::get_as_radian<1>(p3); + CT const lon4 = bg::get_as_radian<0>(p4); + CT const lat4 = bg::get_as_radian<1>(p4); + + CT const distance_result1 = bg::math::hav(lat2 - lat1) + + cos(lat1) * cos(lat2) * bg::math::hav(lon2 - lon1); + CT const distance_result2 = bg::math::hav(lat4 - lat3) + + cos(lat3) * cos(lat4) * bg::math::hav(lon4 - lon3); + + CT const sub = distance_result1 - distance_result2; + + if (sub < -BOOST_GEOMETRY_EPS) + { + result = 1; + } + else if (sub > BOOST_GEOMETRY_EPS) + { + result = 2; + } + else if (fabs(sub) < BOOST_GEOMETRY_EPS) + { + result = 3; + } + return result; + } +}; + +}}} // namespace boost::geometry::formula + +#endif //BOOST_GEOMETRY_FORMULAS_COMPARE_DISTANCE_HAVERSINE_HPP diff --git a/include/boost/geometry/formulas/flat_approximation.hpp b/include/boost/geometry/formulas/flat_approximation.hpp new file mode 100644 index 0000000000..f5891a150d --- /dev/null +++ b/include/boost/geometry/formulas/flat_approximation.hpp @@ -0,0 +1,106 @@ +// Boost.Geometry + +// Copyright (c) 2015-2017 Oracle and/or its affiliates. + +// Contributed and/or modified by Ruoyun Jing, on behalf of Oracle + +// 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_FORMULAS_FLAT_APPROXIMATION_HPP +#define BOOST_GEOMETRY_FORMULAS_FLAT_APPROXIMATION_HPP + +#include +#include +#include + +#include +#include + +#include + +#define BOOST_GEOMETRY_EPS 1e-9 + +namespace bg = boost::geometry; +namespace boost { namespace geometry { namespace formula +{ + +template +class flat_approximation +{ + +public: + typedef int result_type; + template + < + typename T, + typename Spheriod + > + static inline result_type apply(T const& p1, + T const& p2, + T const& p4, + Spheriod const& spheriod) + { + result_type result; + + CT const lon1 = bg::get_as_radian<0>(p1); + CT const lat1 = bg::get_as_radian<1>(p1); + CT const lon2 = bg::get_as_radian<0>(p2); + CT const lat2 = bg::get_as_radian<1>(p2); + CT const lon4 = bg::get_as_radian<0>(p4); + CT const lat4 = bg::get_as_radian<1>(p4); + + CT const distance_result1 = flat_distance(lon1, lat1, lon2, lat2, spheriod); + CT const distance_result2 = flat_distance(lon1, lat1, lon4, lat4, spheriod); + + CT const sub = distance_result1 - distance_result2; + + if (sub < -BOOST_GEOMETRY_EPS) + { + result = 1; + } + else if (sub > BOOST_GEOMETRY_EPS) + { + result = 2; + } + else if (fabs(sub) < BOOST_GEOMETRY_EPS) + { + result = 3; + } + return result; + } +private: + static inline CT flat_distance(CT const& lo1, + CT const& la1, + CT const& lo2, + CT const& la2, + Spheriod const& spher) + { + CT const c2 = CT(2); + CT const earth_r = spher.get_radius(); + CT const earth_f = formula::flattening(spher); + CT const earth_e2 = earth_f * (c2 - earth_f); + + CT const dlat = la2 - la1; + CT const dlon = lo2 - lo1; + + CT const R1 = earth_r * (1 - earth_e2) / + bg::math::sqrt((1 - earth_e2 * bg::math::sqr(sin(la1))) + * (1 - earth_e2 * bg::math::sqr(sin(la1))) + * (1 - earth_e2 * bg::math::sqr(sin(la1)))); + CT const R2 = earth_r / bg::math::sqrt(1 - earth_e2 * bg::math::sqr(sin(la1))); + + CT const dis_North = R1 * dlat; + CT const dis_East = R2 * cos(la1) * dlon; + + CT const distance_result = bg::math::sqrt(bg::math::sqr(dis_North) + + bg::math::sqr(dis_East)); + return distance_result; + } + +}; + +}}} // namespace boost::geometry::formula + +#endif //BOOST_GEOMETRY_FORMULAS_FLAT_APPROXIMATION_HPP