123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266 |
- // Boost.Geometry (aka GGL, Generic Geometry Library)
- // Unit Test
- // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
- // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
- // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
- // 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)
- #include <geometry_test_common.hpp>
- #include <boost/concept/requires.hpp>
- #include <boost/concept_check.hpp>
- #include <boost/core/ignore_unused.hpp>
- #include <boost/geometry/algorithms/assign.hpp>
- #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
- #include <boost/geometry/strategies/concepts/distance_concept.hpp>
- #include <boost/geometry/geometries/point.hpp>
- #ifdef HAVE_TTMATH
- # include <boost/geometry/extensions/contrib/ttmath_stub.hpp>
- #endif
- double const average_earth_radius = 6372795.0;
- template <typename Point, typename LatitudePolicy>
- struct test_distance
- {
- typedef bg::strategy::distance::haversine<double> haversine_type;
- typedef typename bg::strategy::distance::services::return_type<haversine_type, Point, Point>::type return_type;
- BOOST_CONCEPT_ASSERT
- (
- (bg::concepts::PointDistanceStrategy<haversine_type, Point, Point>)
- );
- static void test(double lon1, double lat1, double lon2, double lat2,
- double radius, double expected, double tolerance)
- {
- haversine_type strategy(radius);
- Point p1, p2;
- bg::assign_values(p1, lon1, LatitudePolicy::apply(lat1));
- bg::assign_values(p2, lon2, LatitudePolicy::apply(lat2));
- return_type d = strategy.apply(p1, p2);
- BOOST_CHECK_CLOSE(d, expected, tolerance);
- }
- };
- template <typename Point, typename LatitudePolicy>
- void test_all()
- {
- // earth to unit-sphere -> divide by earth circumference, then it is from 0-1,
- // then multiply with 2 PI, so effectively just divide by earth radius
- double e2u = 1.0 / average_earth_radius;
- // ~ Amsterdam/Paris, 467 kilometers
- double const a_p = 467.2704 * 1000.0;
- test_distance<Point, LatitudePolicy>::test(4, 52, 2, 48, average_earth_radius, a_p, 1.0);
- test_distance<Point, LatitudePolicy>::test(2, 48, 4, 52, average_earth_radius, a_p, 1.0);
- test_distance<Point, LatitudePolicy>::test(4, 52, 2, 48, 1.0, a_p * e2u, 0.001);
- // ~ Amsterdam/Barcelona
- double const a_b = 1232.9065 * 1000.0;
- test_distance<Point, LatitudePolicy>::test(4, 52, 2, 41, average_earth_radius, a_b, 1.0);
- test_distance<Point, LatitudePolicy>::test(2, 41, 4, 52, average_earth_radius, a_b, 1.0);
- test_distance<Point, LatitudePolicy>::test(4, 52, 2, 41, 1.0, a_b * e2u, 0.001);
- }
- template <typename P1, typename P2, typename CalculationType, typename LatitudePolicy>
- void test_services()
- {
- namespace bgsd = bg::strategy::distance;
- namespace services = bg::strategy::distance::services;
- {
- // Compile-check if there is a strategy for this type
- typedef typename services::default_strategy
- <
- bg::point_tag, bg::point_tag, P1, P2
- >::type haversine_strategy_type;
- boost::ignore_unused<haversine_strategy_type>();
- }
- P1 p1;
- bg::assign_values(p1, 4, 52);
- P2 p2;
- bg::assign_values(p2, 2, 48);
- // ~ Amsterdam/Paris, 467 kilometers
- double const km = 1000.0;
- double const a_p = 467.2704 * km;
- double const expected = a_p;
- double const expected_lower = 460.0 * km;
- double const expected_higher = 470.0 * km;
- // 1: normal, calculate distance:
- typedef bgsd::haversine<double, CalculationType> strategy_type;
- typedef typename bgsd::services::return_type<strategy_type, P1, P2>::type return_type;
- strategy_type strategy(average_earth_radius);
- return_type result = strategy.apply(p1, p2);
- BOOST_CHECK_CLOSE(result, return_type(expected), 0.001);
- // 2: the strategy should return the same result if we reverse parameters
- result = strategy.apply(p2, p1);
- BOOST_CHECK_CLOSE(result, return_type(expected), 0.001);
- // 3: "comparable" to construct a "comparable strategy" for P1/P2
- // a "comparable strategy" is a strategy which does not calculate the exact distance, but
- // which returns results which can be mutually compared (e.g. avoid sqrt)
- // 3a: "comparable_type"
- typedef typename services::comparable_type<strategy_type>::type comparable_type;
- // 3b: "get_comparable"
- comparable_type comparable = bgsd::services::get_comparable<strategy_type>::apply(strategy);
- // Check vice versa:
- // First the result of the comparable strategy
- return_type c_result = comparable.apply(p1, p2);
- // Second the comparable result of the expected distance
- return_type c_expected = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected);
- // And that one should be equa.
- BOOST_CHECK_CLOSE(c_result, return_type(c_expected), 0.001);
- // 4: the comparable_type should have a distance_strategy_constructor as well,
- // knowing how to compare something with a fixed distance
- return_type c_dist_lower = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected_lower);
- return_type c_dist_higher = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected_higher);
- // If this is the case:
- BOOST_CHECK(c_dist_lower < c_result && c_result < c_dist_higher);
- // Calculate the Haversine by hand here:
- return_type c_check = return_type(2.0) * asin(sqrt(c_result)) * average_earth_radius;
- BOOST_CHECK_CLOSE(c_check, expected, 0.001);
- // This should also be the case
- return_type dist_lower = services::result_from_distance<strategy_type, P1, P2>::apply(strategy, expected_lower);
- return_type dist_higher = services::result_from_distance<strategy_type, P1, P2>::apply(strategy, expected_higher);
- BOOST_CHECK(dist_lower < result && result < dist_higher);
- }
- /****
- template <typename P, typename Strategy>
- void time_compare_s(int const n)
- {
- boost::timer t;
- P p1, p2;
- bg::assign_values(p1, 1, 1);
- bg::assign_values(p2, 2, 2);
- Strategy strategy;
- typename Strategy::return_type s = 0;
- for (int i = 0; i < n; i++)
- {
- for (int j = 0; j < n; j++)
- {
- s += strategy.apply(p1, p2);
- }
- }
- std::cout << "s: " << s << " t: " << t.elapsed() << std::endl;
- }
- template <typename P>
- void time_compare(int const n)
- {
- time_compare_s<P, bg::strategy::distance::haversine<double> >(n);
- time_compare_s<P, bg::strategy::distance::comparable::haversine<double> >(n);
- }
- #include <time.h>
- double time_sqrt(int n)
- {
- clock_t start = clock();
- double v = 2.0;
- double s = 0.0;
- for (int i = 0; i < n; i++)
- {
- for (int j = 0; j < n; j++)
- {
- s += sqrt(v);
- v += 1.0e-10;
- }
- }
- clock_t end = clock();
- double diff = double(end - start) / CLOCKS_PER_SEC;
- std::cout << "Check: " << s << " Time: " << diff << std::endl;
- return diff;
- }
- double time_normal(int n)
- {
- clock_t start = clock();
- double v = 2.0;
- double s = 0.0;
- for (int i = 0; i < n; i++)
- {
- for (int j = 0; j < n; j++)
- {
- s += v;
- v += 1.0e-10;
- }
- }
- clock_t end = clock();
- double diff = double(end - start) / CLOCKS_PER_SEC;
- std::cout << "Check: " << s << " Time: " << diff << std::endl;
- return diff;
- }
- ***/
- int test_main(int, char* [])
- {
- test_all<bg::model::point<int, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
- test_all<bg::model::point<float, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
- test_all<bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
- // NYI: haversine for mathematical spherical coordinate systems
- // test_all<bg::model::point<double, 2, bg::cs::spherical<bg::degree> >, mathematical_policy>();
- //double t1 = time_sqrt(20000);
- //double t2 = time_normal(20000);
- //std::cout << "Factor: " << (t1 / t2) << std::endl;
- //time_compare<bg::model::point<double, 2, bg::cs::spherical<bg::radian> > >(10000);
- #if defined(HAVE_TTMATH)
- typedef ttmath::Big<1,4> tt;
- test_all<bg::model::point<tt, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
- #endif
- test_services
- <
- bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >,
- bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >,
- double,
- geographic_policy
- >();
- return 0;
- }
|