From c76cee89328c06778b8a27cdf845f8dd3f8dfa79 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Sun, 5 Jun 2011 14:43:07 +0000 Subject: [PATCH] Boost.Geometry: Merged r72086 through r72406 [SVN r72409] --- doc/src/examples/algorithms/area.cpp | 4 +- .../algorithms/area_with_strategy.cpp | 4 +- .../algorithms/length_with_strategy.cpp | 2 +- doc/src/examples/algorithms/transform.cpp | 9 +- doc/src/examples/core/degree_radian.cpp | 4 +- doc/src/examples/quick_start.cpp | 2 +- example/06_a_transformation_example.cpp | 2 +- example/07_a_graph_route_example.cpp | 2 +- example/07_b_graph_route_example.cpp | 4 +- include/boost/geometry/core/cs.hpp | 30 +++- include/boost/geometry/core/tags.hpp | 8 +- .../agnostic/hull_graham_andrew.hpp | 4 +- .../point_in_poly_oriented_winding.hpp | 2 +- .../agnostic/point_in_poly_winding.hpp | 10 +- .../cartesian/distance_projected_point.hpp | 25 ++- .../strategies/cartesian/side_by_triangle.hpp | 21 ++- .../geometry/strategies/intersection.hpp | 2 +- include/boost/geometry/strategies/side.hpp | 19 ++- .../strategies/spherical/area_huiller.hpp | 68 ++++++--- .../strategies/spherical/compare_circular.hpp | 4 +- .../spherical/distance_cross_track.hpp | 61 +++++++- .../spherical/distance_haversine.hpp | 7 +- .../spherical/side_by_cross_track.hpp | 17 +-- .../geometry/strategies/spherical/ssf.hpp | 137 +++++++++++++++++ .../boost/geometry/strategies/strategies.hpp | 1 + .../strategies/strategy_transform.hpp | 97 +++++++++--- test/algorithms/area.cpp | 102 ++++++++++++- test/algorithms/overlay/ccw_traverse.cpp | 2 +- test/algorithms/overlay/traverse.cpp | 2 +- test/algorithms/overlay/traverse_gmp.cpp | 2 +- test/algorithms/simplify.cpp | 4 +- test/algorithms/test_area.hpp | 1 + test/algorithms/within.cpp | 48 ++++++ test/geometry_test_common.hpp | 23 +++ test/strategies/Jamfile.v2 | 1 + test/strategies/cross_track.cpp | 26 ++-- test/strategies/haversine.cpp | 40 ++--- test/strategies/segment_intersection.cpp | 2 +- .../segment_intersection_collinear.cpp | 2 +- test/strategies/side_by_cross_track.cpp | 54 ------- test/strategies/spherical_side.cpp | 142 ++++++++++++++++++ ...oss_track.vcproj => spherical_side.vcproj} | 10 +- test/strategies/strategies_tests.sln | 2 +- 43 files changed, 790 insertions(+), 219 deletions(-) create mode 100644 include/boost/geometry/strategies/spherical/ssf.hpp delete mode 100644 test/strategies/side_by_cross_track.cpp create mode 100644 test/strategies/spherical_side.cpp rename test/strategies/{side_by_cross_track.vcproj => spherical_side.vcproj} (92%) diff --git a/doc/src/examples/algorithms/area.cpp b/doc/src/examples/algorithms/area.cpp index 5b5ef9179..e9dfbbf33 100644 --- a/doc/src/examples/algorithms/area.cpp +++ b/doc/src/examples/algorithms/area.cpp @@ -27,8 +27,8 @@ int main() double area = bg::area(poly); std::cout << "Area: " << area << std::endl; - // Calculate the area of a spherical polygon - bg::model::polygon > > sph_poly; + // Calculate the area of a spherical equatorial polygon + bg::model::polygon > > sph_poly; bg::read_wkt("POLYGON((0 0,0 45,45 0,0 0))", sph_poly); area = bg::area(sph_poly); std::cout << "Area: " << area << std::endl; diff --git a/doc/src/examples/algorithms/area_with_strategy.cpp b/doc/src/examples/algorithms/area_with_strategy.cpp index d47120889..027889a15 100644 --- a/doc/src/examples/algorithms/area_with_strategy.cpp +++ b/doc/src/examples/algorithms/area_with_strategy.cpp @@ -27,8 +27,8 @@ int main() double area = bg::area(poly); std::cout << "Area: " << area << std::endl; - // Calculate the area of a spherical polygon - bg::model::polygon > > sph_poly; + // Calculate the area of a spherical polygon (for latitude: 0 at equator) + bg::model::polygon > > sph_poly; bg::read_wkt("POLYGON((0 0,0 45,45 0,0 0))", sph_poly); area = bg::area(sph_poly); std::cout << "Area: " << area << std::endl; diff --git a/doc/src/examples/algorithms/length_with_strategy.cpp b/doc/src/examples/algorithms/length_with_strategy.cpp index 2c68f3ae8..564d6a846 100644 --- a/doc/src/examples/algorithms/length_with_strategy.cpp +++ b/doc/src/examples/algorithms/length_with_strategy.cpp @@ -17,7 +17,7 @@ int main() { using namespace boost::geometry; - typedef model::point > P; + typedef model::point > P; model::linestring

line; line.push_back(P(2, 41)); line.push_back(P(2, 48)); diff --git a/doc/src/examples/algorithms/transform.cpp b/doc/src/examples/algorithms/transform.cpp index feee4ab13..b2fccaad0 100644 --- a/doc/src/examples/algorithms/transform.cpp +++ b/doc/src/examples/algorithms/transform.cpp @@ -18,7 +18,8 @@ int main() { namespace bg = boost::geometry; - bg::model::point > p1(5.0, 52.0); + // Select a point near the pole (theta=5.0, phi=15.0) + bg::model::point > p1(15.0, 5.0); // Transform from degree to radian. Default strategy is automatically selected, // it will convert from degree to radian @@ -45,9 +46,9 @@ int main() /*` Output: [pre -p1: (5, 52) -p2: (0.0872665, 0.907571) -p3: (0.785012, 0.0686797, 0.615661) +p1: (15, 5) +p2: (0.261799, 0.0872665) +p3: (0.084186, 0.0225576, 0.996195) ] */ //] diff --git a/doc/src/examples/core/degree_radian.cpp b/doc/src/examples/core/degree_radian.cpp index 3cd2d3dba..af27a8652 100644 --- a/doc/src/examples/core/degree_radian.cpp +++ b/doc/src/examples/core/degree_radian.cpp @@ -17,8 +17,8 @@ using namespace boost::geometry; int main() { - typedef model::point > degree_point; - typedef model::point > radian_point; + typedef model::point > degree_point; + typedef model::point > radian_point; degree_point d(4.893, 52.373); radian_point r(0.041, 0.8527); diff --git a/doc/src/examples/quick_start.cpp b/doc/src/examples/quick_start.cpp index ff96b19c7..7f16a7653 100644 --- a/doc/src/examples/quick_start.cpp +++ b/doc/src/examples/quick_start.cpp @@ -135,7 +135,7 @@ int main(void) //[quick_start_spherical typedef boost::geometry::model::point < - double, 2, boost::geometry::cs::spherical + double, 2, boost::geometry::cs::spherical_equatorial > spherical_point; spherical_point amsterdam(4.90, 52.37); diff --git a/example/06_a_transformation_example.cpp b/example/06_a_transformation_example.cpp index 75440f42e..a0f0d2d71 100644 --- a/example/06_a_transformation_example.cpp +++ b/example/06_a_transformation_example.cpp @@ -50,7 +50,7 @@ int main() // - from Cartesian to Spherical coordinate systems and back // - from Cartesian to Cartesian (mapping, affine transformations) and back (inverse) // - Map Projections - // - from Degree to Radian and back in spherical or geographic coordinate systems + // - from Degree to Radian and back in spherical_equatorial or geographic coordinate systems return 0; } diff --git a/example/07_a_graph_route_example.cpp b/example/07_a_graph_route_example.cpp index 929af00a1..8ad672431 100644 --- a/example/07_a_graph_route_example.cpp +++ b/example/07_a_graph_route_example.cpp @@ -267,7 +267,7 @@ int main() // (geographic calculations are in an extension; for sample it makes no difference) typedef boost::geometry::model::point < - double, 2, boost::geometry::cs::spherical + double, 2, boost::geometry::cs::spherical_equatorial > point_type; typedef boost::geometry::model::linestring line_type; diff --git a/example/07_b_graph_route_example.cpp b/example/07_b_graph_route_example.cpp index 17dbea3a6..c28b3c06e 100644 --- a/example/07_b_graph_route_example.cpp +++ b/example/07_b_graph_route_example.cpp @@ -249,11 +249,11 @@ inline void build_route(Graph const& graph, int main() { - // Define a point in the Geographic coordinate system (currently Spherical) + // Define a point in the Geographic coordinate system (currently spherical-equatorial) // (geographic calculations are in an extension; for sample it makes no difference) typedef boost::geometry::model::point < - double, 2, boost::geometry::cs::spherical + double, 2, boost::geometry::cs::spherical_equatorial > point_type; typedef boost::geometry::model::linestring line_type; diff --git a/include/boost/geometry/core/cs.hpp b/include/boost/geometry/core/cs.hpp index 93e0eec77..a53d6e278 100644 --- a/include/boost/geometry/core/cs.hpp +++ b/include/boost/geometry/core/cs.hpp @@ -79,7 +79,7 @@ struct geographic /*! -\brief Spherical coordinate system, in degree or in radian +\brief Spherical (polar) coordinate system, in degree or in radian \details Defines the spherical coordinate system where points are defined in two angles and an optional radius usually known as r, theta, phi @@ -102,6 +102,25 @@ struct spherical typedef DegreeOrRadian units; }; + +/*! +\brief Spherical equatorial coordinate system, in degree or in radian +\details This one resembles the geographic coordinate system, and has latitude + up from zero at the equator, to 90 at the pole + (opposite to the spherical(polar) coordinate system). + Used in astronomy and in GIS (but there is also the geographic) + +\see http://en.wikipedia.org/wiki/Spherical_coordinates +\ingroup cs +*/ +template +struct spherical_equatorial +{ + typedef DegreeOrRadian units; +}; + + + /*! \brief Polar coordinate system \details Defines the polar coordinate system "in which each point @@ -144,9 +163,16 @@ struct cs_tag > template struct cs_tag > { - typedef spherical_tag type; + typedef spherical_polar_tag type; }; +template +struct cs_tag > +{ + typedef spherical_equatorial_tag type; +}; + + template<> struct cs_tag { diff --git a/include/boost/geometry/core/tags.hpp b/include/boost/geometry/core/tags.hpp index 6469b6c11..f110a079a 100644 --- a/include/boost/geometry/core/tags.hpp +++ b/include/boost/geometry/core/tags.hpp @@ -24,11 +24,15 @@ namespace boost { namespace geometry /// Tag indicating Cartesian coordinate system family (cartesian,epsg) struct cartesian_tag {}; +/// Tag indicating Spherical polar coordinate system family +struct spherical_polar_tag {}; + +/// Tag indicating Spherical equatorial coordinate system family +struct spherical_equatorial_tag {}; + /// Tag indicating Geographic coordinate system family (geographic) struct geographic_tag {}; -/// Tag indicating Spherical coordinate system family (spherical,celestial,...) -struct spherical_tag {}; // Tags defining tag hierarchy diff --git a/include/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp b/include/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp index 1b91a1c23..037e390b5 100644 --- a/include/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp +++ b/include/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp @@ -301,7 +301,7 @@ public: range_type, range_iterator, container_type, - typename strategy_side::type + typename strategy::side::services::default_strategy::type > assigner(extremes.left, extremes.right); geometry::detail::for_each_range(geometry, assigner); @@ -358,7 +358,7 @@ private: template static inline void add_to_hull(point_type const& p, container_type& output) { - typedef typename strategy_side::type side; + typedef typename strategy::side::services::default_strategy::type side; output.push_back(p); register std::size_t output_size = output.size(); diff --git a/include/boost/geometry/strategies/agnostic/point_in_poly_oriented_winding.hpp b/include/boost/geometry/strategies/agnostic/point_in_poly_oriented_winding.hpp index bb550e61b..a7ebc7f75 100644 --- a/include/boost/geometry/strategies/agnostic/point_in_poly_oriented_winding.hpp +++ b/include/boost/geometry/strategies/agnostic/point_in_poly_oriented_winding.hpp @@ -62,7 +62,7 @@ class oriented_winding >::type calculation_type; - typedef typename strategy_side + typedef typename strategy::side::services::default_strategy < typename cs_tag::type >::type strategy_side_type; diff --git a/include/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp b/include/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp index 373014f15..c98de8ab3 100644 --- a/include/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp +++ b/include/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp @@ -58,7 +58,7 @@ class winding >::type calculation_type; - typedef typename strategy_side + typedef typename strategy::side::services::default_strategy < typename cs_tag::type >::type strategy_side_type; @@ -188,7 +188,13 @@ struct default_strategy -struct default_strategy +struct default_strategy +{ + typedef winding type; +}; + +template +struct default_strategy { typedef winding type; }; diff --git a/include/boost/geometry/strategies/cartesian/distance_projected_point.hpp b/include/boost/geometry/strategies/cartesian/distance_projected_point.hpp index 4378949ea..063571fc3 100644 --- a/include/boost/geometry/strategies/cartesian/distance_projected_point.hpp +++ b/include/boost/geometry/strategies/cartesian/distance_projected_point.hpp @@ -116,13 +116,14 @@ public : { assert_dimension_equal(); - /* Algorithm - POINT v(x2 - x1, y2 - y1); - POINT w(px - x1, py - y1); - c1 = w . v - c2 = v . v - b = c1 / c2 - RETURN POINT(x1 + b * vx, y1 + b * vy); + /* + Algorithm [p1: (x1,y1), p2: (x2,y2), p: (px,py)] + VECTOR v(x2 - x1, y2 - y1) + VECTOR w(px - x1, py - y1) + c1 = w . v + c2 = v . v + b = c1 / c2 + RETURN POINT(x1 + b * vx, y1 + b * vy) */ // v is multiplied below with a (possibly) FP-value, so should be in FP @@ -137,21 +138,20 @@ public : Strategy strategy; boost::ignore_unused_variable_warning(strategy); - calculation_type zero = calculation_type(); - fp_type c1 = dot_product(w, v); + calculation_type const zero = calculation_type(); + fp_type const c1 = dot_product(w, v); if (c1 <= zero) { return strategy.apply(p, p1); } - fp_type c2 = dot_product(v, v); + fp_type const c2 = dot_product(v, v); if (c2 <= c1) { return strategy.apply(p, p2); } // See above, c1 > 0 AND c2 > c1 so: c2 != 0 - fp_type b = fp_type(c1) / fp_type(c2); - + fp_type const b = c1 / c2; fp_strategy_type fp_strategy = strategy::distance::services::get_similar @@ -168,7 +168,6 @@ public : return fp_strategy.apply(p, projected); } - }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS diff --git a/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp b/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp index 123689975..da3c3ed77 100644 --- a/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp +++ b/include/boost/geometry/strategies/cartesian/side_by_triangle.hpp @@ -21,6 +21,7 @@ #include +#include namespace boost { namespace geometry @@ -90,22 +91,28 @@ public : promoted_type const s = dx * dpy - dy * dpx; promoted_type zero = promoted_type(); - return math::equals(s, zero) ? 0 : s > zero ? 1 : -1; - //return s > 0 ? 1 : s < 0 ? -1 : 0; + return math::equals(s, zero) ? 0 + : s > zero ? 1 + : -1; } }; -}} // namespace strategy::side - #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS -template -struct strategy_side +namespace services { - typedef strategy::side::side_by_triangle type; + +template +struct default_strategy +{ + typedef side_by_triangle type; }; + +} #endif +}} // namespace strategy::side + }} // namespace boost::geometry diff --git a/include/boost/geometry/strategies/intersection.hpp b/include/boost/geometry/strategies/intersection.hpp index 056910191..234d845d6 100644 --- a/include/boost/geometry/strategies/intersection.hpp +++ b/include/boost/geometry/strategies/intersection.hpp @@ -77,7 +77,7 @@ public: CalculationType > segment_intersection_strategy_type; - typedef typename strategy_side + typedef typename strategy::side::services::default_strategy < Tag, CalculationType diff --git a/include/boost/geometry/strategies/side.hpp b/include/boost/geometry/strategies/side.hpp index 55a168b77..31587a87c 100644 --- a/include/boost/geometry/strategies/side.hpp +++ b/include/boost/geometry/strategies/side.hpp @@ -21,6 +21,11 @@ namespace boost { namespace geometry { +namespace strategy { namespace side +{ + +namespace services +{ /*! \brief Traits class binding a side determination strategy to a coordinate system @@ -29,12 +34,22 @@ namespace boost { namespace geometry \tparam CalculationType \tparam_calculation */ template -struct strategy_side +struct default_strategy { - typedef strategy::not_implemented type; + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types) + ); }; +} // namespace services + + +}} // namespace strategy::side + + }} // namespace boost::geometry #endif // BOOST_GEOMETRY_STRATEGIES_SIDE_HPP diff --git a/include/boost/geometry/strategies/spherical/area_huiller.hpp b/include/boost/geometry/strategies/spherical/area_huiller.hpp index 23fd0cd64..0b739dcef 100644 --- a/include/boost/geometry/strategies/spherical/area_huiller.hpp +++ b/include/boost/geometry/strategies/spherical/area_huiller.hpp @@ -10,7 +10,6 @@ #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP -#include #include @@ -45,6 +44,8 @@ http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/40409/1/07-03.pdf, is simple and works well in most cases but not in 180 meridian crossing cases. This probably could be solved. +\note This version is made for spherical equatorial coordinate systems + \qbk{ [heading Example] @@ -64,10 +65,21 @@ template > class huiller { +typedef typename boost::mpl::if_c + < + boost::is_void::type::value, + typename select_most_precise + < + typename coordinate_type::type, + double + >::type, + CalculationType + >::type calculation_type; + protected : struct excess_sum { - double sum; + calculation_type sum; // Distances are calculated on unit sphere here strategy::distance::haversine @@ -78,18 +90,18 @@ protected : : sum(0) , distance_over_unit_sphere(1) {} - inline double area(double radius) const + inline calculation_type area(calculation_type radius) const { return - sum * radius * radius; } }; public : - typedef double return_type; + typedef calculation_type return_type; typedef PointOfSegment segment_point_type; typedef excess_sum state_type; - inline huiller(double radius = 1.0) + inline huiller(calculation_type radius = 1.0) : m_radius(radius) {} @@ -99,26 +111,25 @@ public : { if (! geometry::math::equals(get<0>(p1), get<0>(p2))) { - namespace mc = boost::math::constants; - - double const two_pi = 2.0 * mc::pi(); - double const half = 0.5; - double const two = 2.0; - double const four = 4.0; + calculation_type const half = 0.5; + calculation_type const two = 2.0; + calculation_type const four = 4.0; + calculation_type const two_pi = two * geometry::math::pi(); + calculation_type const half_pi = half * geometry::math::pi(); // Distance p1 p2 - double a = state.distance_over_unit_sphere.apply(p1, p2); + calculation_type a = state.distance_over_unit_sphere.apply(p1, p2); // Sides on unit sphere to south pole - double b = half * mc::pi() - geometry::get_as_radian<1>(p2); - double c = half * mc::pi() - geometry::get_as_radian<1>(p1); + calculation_type b = half_pi - geometry::get_as_radian<1>(p2); + calculation_type c = half_pi - geometry::get_as_radian<1>(p1); // Semi parameter - double s = half * (a + b + c); + calculation_type s = half * (a + b + c); // E: spherical excess, using l'Huiller's formula // [tg(e / 4)]2 = tg[s / 2] tg[(s-a) / 2] tg[(s-b) / 2] tg[(s-c) / 2] - double E = four * atan(sqrt(geometry::math::abs(tan(s / two) + calculation_type E = four * atan(sqrt(geometry::math::abs(tan(s / two) * tan((s - a) / two) * tan((s - b) / two) * tan((s - c) / two)))); @@ -130,11 +141,11 @@ public : // we have to take the dateline into account. // TODO: check this / enhance this, should be more robust. See also the "grow" for ll // TODO: use minmax or "smaller"/"compare" strategy for this - double lon1 = geometry::get_as_radian<0>(p1) < 0 + calculation_type lon1 = geometry::get_as_radian<0>(p1) < 0 ? geometry::get_as_radian<0>(p1) + two_pi : geometry::get_as_radian<0>(p1); - double lon2 = geometry::get_as_radian<0>(p2) < 0 + calculation_type lon2 = geometry::get_as_radian<0>(p2) < 0 ? geometry::get_as_radian<0>(p2) + two_pi : geometry::get_as_radian<0>(p2); @@ -154,18 +165,27 @@ public : private : /// Radius of the sphere - double m_radius; + calculation_type m_radius; }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace services { - template - struct default_strategy - { - typedef strategy::area::huiller type; - }; + + +template +struct default_strategy +{ + typedef strategy::area::huiller type; +}; + +// Note: spherical polar coordinate system requires "get_as_radian_equatorial" +/***template +struct default_strategy +{ + typedef strategy::area::huiller type; +};***/ } // namespace services diff --git a/include/boost/geometry/strategies/spherical/compare_circular.hpp b/include/boost/geometry/strategies/spherical/compare_circular.hpp index e33147c22..fee1e2b7e 100644 --- a/include/boost/geometry/strategies/spherical/compare_circular.hpp +++ b/include/boost/geometry/strategies/spherical/compare_circular.hpp @@ -117,7 +117,7 @@ template template class CoordinateSystem, typename Units > -struct strategy_compare, 0> +struct strategy_compare, 0> { typedef typename coordinate_type::type coordinate_type; typedef strategy::compare::circular_comparator @@ -134,7 +134,7 @@ template template class CoordinateSystem, typename Units > -struct strategy_compare, 0> +struct strategy_compare, 0> { typedef typename coordinate_type::type coordinate_type; typedef strategy::compare::circular_comparator diff --git a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp index 5c5f3e9a2..f902591bf 100644 --- a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -219,7 +219,12 @@ struct comparable_type +template +< + typename Point, typename PointOfSegment, + typename CalculationType, + typename Strategy +> struct get_comparable > { typedef typename comparable_type @@ -234,7 +239,12 @@ public : }; -template +template +< + typename Point, typename PointOfSegment, + typename CalculationType, + typename Strategy +> struct result_from_distance > { private : @@ -248,7 +258,12 @@ public : }; -template +template +< + typename Point, typename PointOfSegment, + typename CalculationType, + typename Strategy +> struct strategy_point_point > { typedef Strategy type; @@ -256,9 +271,17 @@ struct strategy_point_point" template -struct default_strategy +struct default_strategy + < + segment_tag, Point, PointOfSegment, + spherical_polar_tag, spherical_polar_tag, + Strategy + > { typedef cross_track < @@ -271,12 +294,40 @@ struct default_strategy::type, Strategy >::type > type; }; +*/ + +template +struct default_strategy + < + segment_tag, Point, PointOfSegment, + spherical_equatorial_tag, spherical_equatorial_tag, + Strategy + > +{ + typedef cross_track + < + Point, + PointOfSegment, + void, + typename boost::mpl::if_ + < + boost::is_void, + typename default_strategy + < + point_tag, Point, PointOfSegment, + spherical_equatorial_tag, spherical_equatorial_tag + >::type, + Strategy + >::type + > type; +}; + } // namespace services diff --git a/include/boost/geometry/strategies/spherical/distance_haversine.hpp b/include/boost/geometry/strategies/spherical/distance_haversine.hpp index 3722de23c..eb9bf553b 100644 --- a/include/boost/geometry/strategies/spherical/distance_haversine.hpp +++ b/include/boost/geometry/strategies/spherical/distance_haversine.hpp @@ -306,13 +306,16 @@ public : }; -// Register it as the default for point-types in a spherical coordinate system +// Register it as the default for point-types +// in a spherical equatorial coordinate system template -struct default_strategy +struct default_strategy { typedef strategy::distance::haversine type; }; +// Note: spherical polar coordinate system requires "get_as_radian_equatorial" + } // namespace services #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS diff --git a/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp b/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp index 4342fc6dd..bbae9df90 100644 --- a/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp +++ b/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp @@ -83,31 +83,18 @@ public : CalculationType >::type coordinate_type; - // Calculate the distance using the Haversine formula. - // That is also applicable on the spherical earth. A radius is not necessary. - double d1 = 0.001; // m_strategy.apply(sp1, p); double crs_AD = detail::course(p1, p); double crs_AB = detail::course(p1, p2); - double XTD = geometry::math::abs(asin(sin(d1) * sin(crs_AD - crs_AB))); + double XTD = asin(sin(d1) * sin(crs_AD - crs_AB)); - return math::equals(XTD, 0) ? 0 : XTD > 0 ? 1 : -1; - //return s > 0 ? 1 : s < 0 ? -1 : 0; + return math::equals(XTD, 0) ? 0 : XTD < 0 ? 1 : -1; } }; }} // namespace strategy::side -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS -template -struct strategy_side -{ - typedef strategy::side::side_by_cross_track type; -}; -#endif - - }} // namespace boost::geometry diff --git a/include/boost/geometry/strategies/spherical/ssf.hpp b/include/boost/geometry/strategies/spherical/ssf.hpp new file mode 100644 index 000000000..96a7c379d --- /dev/null +++ b/include/boost/geometry/strategies/spherical/ssf.hpp @@ -0,0 +1,137 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2011 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_SPHERICAL_SSF_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP + +#include +#include + +#include +#include +#include + +#include +#include + +#include +//#include + + +namespace boost { namespace geometry +{ + + +namespace strategy { namespace side +{ + + +/*! +\brief Check at which side of a segment a point lies: +\details from a Great Circle segment between two points: + left of segment (> 0), right of segment (< 0), on segment (0) +\ingroup strategies +\tparam CalculationType \tparam_calculation + */ +template +class spherical_side_formula +{ + +public : + template + static inline int apply(P1 const& p1, P2 const& p2, P const& p) + { + typedef typename boost::mpl::if_c + < + boost::is_void::type::value, + + // Select at least a double... + typename select_most_precise + < + typename select_most_precise + < + typename select_most_precise + < + typename coordinate_type::type, + typename coordinate_type::type + >::type, + typename coordinate_type

::type + >::type, + double + >::type, + CalculationType + >::type coordinate_type; + + // Convenient shortcuts + typedef coordinate_type ct; + ct const lambda1 = get_as_radian<0>(p1); + ct const delta1 = get_as_radian<1>(p1); + ct const lambda2 = get_as_radian<0>(p2); + ct const delta2 = get_as_radian<1>(p2); + ct const lambda = get_as_radian<0>(p); + ct const delta = get_as_radian<1>(p); + + // Create temporary points (vectors) on unit a sphere + ct const cos_delta1 = cos(delta1); + ct const c1x = cos_delta1 * cos(lambda1); + ct const c1y = cos_delta1 * sin(lambda1); + ct const c1z = sin(delta1); + + ct const cos_delta2 = cos(delta2); + ct const c2x = cos_delta2 * cos(lambda2); + ct const c2y = cos_delta2 * sin(lambda2); + ct const c2z = sin(delta2); + + // (Third point is converted directly) + ct const cos_delta = cos(delta); + + // Apply the "Spherical Side Formula" as presented on my blog + ct const dist + = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda) + + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda) + + (c1x * c2y - c1y * c2x) * sin(delta); + + ct zero = ct(); + return dist > zero ? 1 + : dist < zero ? -1 + : 0; + } +}; + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +/*template +struct default_strategy +{ + typedef spherical_side_formula type; +};*/ + +template +struct default_strategy +{ + typedef spherical_side_formula type; +}; + +template +struct default_strategy +{ + typedef spherical_side_formula type; +}; + +} +#endif + +}} // namespace strategy::side + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SSF_HPP diff --git a/include/boost/geometry/strategies/strategies.hpp b/include/boost/geometry/strategies/strategies.hpp index 30668342d..25369a7bf 100644 --- a/include/boost/geometry/strategies/strategies.hpp +++ b/include/boost/geometry/strategies/strategies.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include diff --git a/include/boost/geometry/strategies/strategy_transform.hpp b/include/boost/geometry/strategies/strategy_transform.hpp index 35c6a4526..34e19fc77 100644 --- a/include/boost/geometry/strategies/strategy_transform.hpp +++ b/include/boost/geometry/strategies/strategy_transform.hpp @@ -155,26 +155,48 @@ namespace detail /// Helper function for conversion, phi/theta are in radians template - inline void spherical_to_cartesian(T phi, T theta, R r, P& p) + inline void spherical_polar_to_cartesian(T phi, T theta, R r, P& p) { assert_dimension(); // http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_spherical_coordinates + // http://www.vias.org/comp_geometry/math_coord_convert_3d.htm + // https://moodle.polymtl.ca/file.php/1183/Autres_Documents/Derivation_for_Spherical_Co-ordinates.pdf + // http://en.citizendium.org/wiki/Spherical_polar_coordinates + // Phi = first, theta is second, r is third, see documentation on cs::spherical // (calculations are splitted to implement ttmath) T r_sin_theta = r; + T r_cos_theta = r; r_sin_theta *= sin(theta); + r_cos_theta *= cos(theta); set<0>(p, r_sin_theta * cos(phi)); set<1>(p, r_sin_theta * sin(phi)); - - T r_cos_theta = r; - r_cos_theta *= cos(theta); - set<2>(p, r_cos_theta); } + + /// Helper function for conversion, lambda/delta (lon lat) are in radians + template + inline void spherical_equatorial_to_cartesian(T lambda, T delta, R r, P& p) + { + assert_dimension(); + + // http://mathworld.wolfram.com/GreatCircle.html + // http://www.spenvis.oma.be/help/background/coortran/coortran.html WRONG + + T r_cos_delta = r; + T r_sin_delta = r; + r_cos_delta *= cos(delta); + r_sin_delta *= sin(delta); + + set<0>(p, r_cos_delta * cos(lambda)); + set<1>(p, r_cos_delta * sin(lambda)); + set<2>(p, r_sin_delta); + } + /// Helper function for conversion template @@ -230,16 +252,28 @@ namespace detail \tparam P2 second point type */ template -struct from_spherical_2_to_cartesian_3 +struct from_spherical_polar_2_to_cartesian_3 { inline bool apply(P1 const& p1, P2& p2) const { assert_dimension(); - detail::spherical_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2); + detail::spherical_polar_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2); return true; } }; +template +struct from_spherical_equatorial_2_to_cartesian_3 +{ + inline bool apply(P1 const& p1, P2& p2) const + { + assert_dimension(); + detail::spherical_equatorial_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2); + return true; + } +}; + + /*! \brief Transformation strategy for 3D spherical (phi,theta,r) to 3D cartesian (x,y,z) \ingroup transform @@ -247,17 +281,30 @@ struct from_spherical_2_to_cartesian_3 \tparam P2 second point type */ template -struct from_spherical_3_to_cartesian_3 +struct from_spherical_polar_3_to_cartesian_3 { inline bool apply(P1 const& p1, P2& p2) const { assert_dimension(); - detail::spherical_to_cartesian( + detail::spherical_polar_to_cartesian( get_as_radian<0>(p1), get_as_radian<1>(p1), get<2>(p1), p2); return true; } }; +template +struct from_spherical_equatorial_3_to_cartesian_3 +{ + inline bool apply(P1 const& p1, P2& p2) const + { + assert_dimension(); + detail::spherical_equatorial_to_cartesian( + get_as_radian<0>(p1), get_as_radian<1>(p1), get<2>(p1), p2); + return true; + } +}; + + /*! \brief Transformation strategy for 3D cartesian (x,y,z) to 2D spherical (phi,theta) \details on Unit sphere @@ -267,7 +314,7 @@ struct from_spherical_3_to_cartesian_3 \note If x,y,z point is not lying on unit sphere, transformation will return false */ template -struct from_cartesian_3_to_spherical_2 +struct from_cartesian_3_to_spherical_polar_2 { inline bool apply(P1 const& p1, P2& p2) const { @@ -284,7 +331,7 @@ struct from_cartesian_3_to_spherical_2 \tparam P2 second point type */ template -struct from_cartesian_3_to_spherical_3 +struct from_cartesian_3_to_spherical_polar_3 { inline bool apply(P1 const& p1, P2& p2) const { @@ -343,30 +390,42 @@ struct default_strategy, CoordSys -struct default_strategy +struct default_strategy { - typedef from_spherical_2_to_cartesian_3 type; + typedef from_spherical_polar_2_to_cartesian_3 type; }; /// Specialization to transform from sphere(phi,theta,r) to XYZ template -struct default_strategy +struct default_strategy { - typedef from_spherical_3_to_cartesian_3 type; + typedef from_spherical_polar_3_to_cartesian_3 type; +}; + +template +struct default_strategy +{ + typedef from_spherical_equatorial_2_to_cartesian_3 type; +}; + +template +struct default_strategy +{ + typedef from_spherical_equatorial_3_to_cartesian_3 type; }; /// Specialization to transform from XYZ to unit sphere(phi,theta) template -struct default_strategy +struct default_strategy { - typedef from_cartesian_3_to_spherical_2 type; + typedef from_cartesian_3_to_spherical_polar_2 type; }; /// Specialization to transform from XYZ to sphere(phi,theta,r) template -struct default_strategy +struct default_strategy { - typedef from_cartesian_3_to_spherical_3 type; + typedef from_cartesian_3_to_spherical_polar_3 type; }; diff --git a/test/algorithms/area.cpp b/test/algorithms/area.cpp index 4ed4e7009..2cf80569f 100644 --- a/test/algorithms/area.cpp +++ b/test/algorithms/area.cpp @@ -13,8 +13,6 @@ // http://www.boost.org/LICENSE_1_0.txt) -#include - #include #include @@ -60,15 +58,19 @@ void test_all() } template -void test_spherical() +void test_spherical(bool polar = false) { + typedef typename bg::coordinate_type::type ct; bg::model::polygon geometry; // unit-sphere has area of 4-PI. Polygon covering 1/8 of it: - double expected = 4.0 * boost::math::constants::pi() / 8.0; + // calculations splitted for ttmath + ct const four = 4.0; + ct const eight = 8.0; + ct expected = four * boost::geometry::math::pi() / eight; bg::read_wkt("POLYGON((0 0,0 90,90 0,0 0))", geometry); - double area = bg::area(geometry); + ct area = bg::area(geometry); BOOST_CHECK_CLOSE(area, expected, 0.0001); // With strategy, radius 2 -> 4 pi r^2 @@ -78,7 +80,91 @@ void test_spherical() > strategy(2.0); area = bg::area(geometry, strategy); - BOOST_CHECK_CLOSE(area, 2.0 * 2.0 * expected, 0.0001); + ct const two = 2.0; + BOOST_CHECK_CLOSE(area, two * two * expected, 0.0001); + + // Wrangel Island (dateline crossing) + // With (spherical) Earth strategy + bg::strategy::area::huiller + < + typename bg::point_type::type + > spherical_earth(6373); + bg::read_wkt("POLYGON((-178.7858 70.7852, 177.4758 71.2333, 179.7436 71.5733, -178.7858 70.7852))", geometry); + area = bg::area(geometry, spherical_earth); + // SQL Server gives: 4537.9654419375 + // PostGIS gives: 4537.9311668307 + // Note: those are Geographic, this test is Spherical + BOOST_CHECK_CLOSE(area, 4506.6389, 0.001); + + // Wrangel, more in detail + bg::read_wkt("POLYGON((-178.568604 71.564148,-178.017548 71.449692,-177.833313 71.3461,-177.502838 71.277466 ,-177.439453 71.226929,-177.620026 71.116638,-177.9389 71.037491,-178.8186 70.979965,-179.274445 70.907761,-180 70.9972,179.678314 70.895538,179.272766 70.888596,178.791016 70.7964,178.617737 71.035538,178.872192 71.217484,179.530273 71.4383 ,-180 71.535843 ,-179.628601 71.577194,-179.305298 71.551361,-179.03421 71.597748,-178.568604 71.564148))", geometry); + area = bg::area(geometry, spherical_earth); + // SQL Server gives: 7669.10402181435 + // PostGIS gives: 7669.55565459832 + BOOST_CHECK_CLOSE(area, 7616.523769, 0.001); + + // Check more at the equator + /* + select 1,geography::STGeomFromText('POLYGON((-178.7858 10.7852 , 179.7436 11.5733 , 177.4758 11.2333 , -178.7858 10.7852))',4326) .STArea()/1000000.0 + union select 2,geography::STGeomFromText('POLYGON((-178.7858 20.7852 , 179.7436 21.5733 , 177.4758 21.2333 , -178.7858 20.7852))',4326) .STArea()/1000000.0 + union select 3,geography::STGeomFromText('POLYGON((-178.7858 30.7852 , 179.7436 31.5733 , 177.4758 31.2333 , -178.7858 30.7852))',4326) .STArea()/1000000.0 + union select 0,geography::STGeomFromText('POLYGON((-178.7858 0.7852 , 179.7436 1.5733 , 177.4758 1.2333 , -178.7858 0.7852))',4326) .STArea()/1000000.0 + union select 4,geography::STGeomFromText('POLYGON((-178.7858 40.7852 , 179.7436 41.5733 , 177.4758 41.2333 , -178.7858 40.7852))',4326) .STArea()/1000000.0 + */ + + bg::read_wkt("POLYGON((-178.7858 0.7852, 177.4758 1.2333, 179.7436 1.5733, -178.7858 0.7852))", geometry); + area = bg::area(geometry, spherical_earth); + BOOST_CHECK_CLOSE(area, 14136.09946, 0.001); // SQL Server gives: 14064.1902284513 + + + bg::read_wkt("POLYGON((-178.7858 10.7852, 177.4758 11.2333, 179.7436 11.5733, -178.7858 10.7852))", geometry); + area = bg::area(geometry, spherical_earth); + BOOST_CHECK_CLOSE(area, 13760.2456, 0.001); // SQL Server gives: 13697.0941155193 + + bg::read_wkt("POLYGON((-178.7858 20.7852, 177.4758 21.2333, 179.7436 21.5733, -178.7858 20.7852))", geometry); + area = bg::area(geometry, spherical_earth); + BOOST_CHECK_CLOSE(area, 12987.8682, 0.001); // SQL Server gives: 12944.3970990317 -> -39m^2 + + bg::read_wkt("POLYGON((-178.7858 30.7852, 177.4758 31.2333, 179.7436 31.5733, -178.7858 30.7852))", geometry); + area = bg::area(geometry, spherical_earth); + BOOST_CHECK_CLOSE(area, 11856.3935, 0.001); // SQL Server gives: 11838.5338423574 -> -18m^2 + + bg::read_wkt("POLYGON((-178.7858 40.7852, 177.4758 41.2333, 179.7436 41.5733, -178.7858 40.7852))", geometry); + area = bg::area(geometry, spherical_earth); + BOOST_CHECK_CLOSE(area, 10404.627685523914, 0.001); // SQL Server gives: 10412.0607137119, -> +8m^2 + + // Concave + bg::read_wkt("POLYGON((0 40,1 42,0 44,2 43,4 44,3 42,4 40,2 41,0 40))", geometry); + area = bg::area(geometry, spherical_earth); + BOOST_CHECK_CLOSE(area, 73538.2958, 0.001); // SQL Server gives: 73604.2047689719 + + // With hole POLYGON((0 40,4 40,4 44,0 44,0 40),(1 41,2 43,3 42,1 41)) + bg::read_wkt("POLYGON((0 40,0 44,4 44,4 40,0 40),(1 41,3 42,2 43,1 41))", geometry); + area = bg::area(geometry, spherical_earth); + BOOST_CHECK_CLOSE(area, 133233.844876, 0.001); // SQL Server gives: 133353.335 + + { + bg::model::ring aurha; // a'dam-utr-rott.-den haag-a'dam + bg::read_wkt("POLYGON((4.892 52.373,5.119 52.093,4.479 51.930,4.23 52.08,4.892 52.373))", aurha); + if (polar) + { + // Create colatitudes (measured from pole) + BOOST_FOREACH(Point& p, aurha) + { + bg::set<1>(p, ct(90) - bg::get<1>(p)); + } + bg::correct(aurha); + } + bg::strategy::area::huiller + < + typename bg::point_type::type + > huiller(6372.795); + area = bg::area(aurha, huiller); + BOOST_CHECK_CLOSE(area, 1476.645675, 0.0001); + + // SQL Server gives: 1481.55595960659 + // for select geography::STGeomFromText('POLYGON((4.892 52.373,4.23 52.08,4.479 51.930,5.119 52.093,4.892 52.373))',4326).STArea()/1000000.0 + } } template @@ -116,7 +202,8 @@ int test_main(int, char* []) test_all >(); test_all >(); - test_spherical > >(); + test_spherical > >(); + //test_spherical > >(true); test_ccw >(); test_open >(); @@ -124,6 +211,7 @@ int test_main(int, char* []) #ifdef HAVE_TTMATH test_all >(); + test_spherical > >(); #endif return 0; diff --git a/test/algorithms/overlay/ccw_traverse.cpp b/test/algorithms/overlay/ccw_traverse.cpp index 82f9d3b1e..463ecdbde 100644 --- a/test/algorithms/overlay/ccw_traverse.cpp +++ b/test/algorithms/overlay/ccw_traverse.cpp @@ -35,7 +35,7 @@ template inline typename bg::coordinate_type::type intersect(Geometry1 const& g1, Geometry2 const& g2, std::string const& name, bg::detail::overlay::operation_type op) { - typedef typename bg::strategy_side + typedef typename bg::strategy::side::services::default_strategy < typename bg::cs_tag::type >::type side_strategy_type; diff --git a/test/algorithms/overlay/traverse.cpp b/test/algorithms/overlay/traverse.cpp index baf0653cb..480b204aa 100644 --- a/test/algorithms/overlay/traverse.cpp +++ b/test/algorithms/overlay/traverse.cpp @@ -134,7 +134,7 @@ struct test_traverse //std::cout << bg::area(g1) << " " << bg::area(g2) << std::endl; #endif - typedef typename bg::strategy_side + typedef typename bg::strategy::side::services::default_strategy < typename bg::cs_tag::type >::type side_strategy_type; diff --git a/test/algorithms/overlay/traverse_gmp.cpp b/test/algorithms/overlay/traverse_gmp.cpp index b8d6a4272..24a4b2116 100644 --- a/test/algorithms/overlay/traverse_gmp.cpp +++ b/test/algorithms/overlay/traverse_gmp.cpp @@ -62,7 +62,7 @@ void test_traverse(std::string const& caseid, G1 const& g1, G2 const& g2) typedef std::vector ip_vector; ip_vector ips; - typedef typename bg::strategy_side + typedef typename bg::strategy::side::services::default_strategy < typename bg::cs_tag::type >::type strategy_type; diff --git a/test/algorithms/simplify.cpp b/test/algorithms/simplify.cpp index 2555a638e..09ca32135 100644 --- a/test/algorithms/simplify.cpp +++ b/test/algorithms/simplify.cpp @@ -96,11 +96,11 @@ int test_main(int, char* []) test_all >(); test_all >(); - test_spherical > >(); + test_spherical > >(); #if defined(HAVE_TTMATH) test_all >(); - test_spherical > >(); + test_spherical > >(); #endif return 0; diff --git a/test/algorithms/test_area.hpp b/test/algorithms/test_area.hpp index 5441bcffc..bf7d7a021 100644 --- a/test/algorithms/test_area.hpp +++ b/test/algorithms/test_area.hpp @@ -14,6 +14,7 @@ #include #include +#include #include #include diff --git a/test/algorithms/within.cpp b/test/algorithms/within.cpp index ad59f4d78..b8820f99b 100644 --- a/test/algorithms/within.cpp +++ b/test/algorithms/within.cpp @@ -62,7 +62,51 @@ void test_all() true, false); } +template +void test_spherical() +{ + typedef typename bg::coordinate_type::type ct; + bg::model::polygon wrangel; + // SQL Server check (no geography::STWithin, so check with intersection trick) + /* + + with q as ( + select geography::STGeomFromText('POLYGON((-178.569 71.5641,-179.034 71.5977,-179.305 71.5514,-179.629 71.5772,-180 71.5358,179.53 71.4383,178.872 71.2175,178.618 71.0355,178.791 70.7964,179.273 70.8886,179.678 70.8955,-180 70.9972,-179.274 70.9078,-178.819 70.98,-177.939 71.0375,-177.62 71.1166,-177.439 71.2269,-177.503 71.2775,-177.833 71.3461,-178.018 71.4497,-178.569 71.5641))',4326) as wrangel + ) + + select wrangel.STArea()/1000000.0 + ,geography::STGeomFromText('POINT(-179.3 71.27)',4326).STIntersection(wrangel).STAsText() as workaround_within_1 + ,geography::STGeomFromText('POINT(-179.9 70.95)',4326).STIntersection(wrangel).STAsText() as workaround_within_2 + ,geography::STGeomFromText('POINT(179.9 70.95)',4326).STIntersection(wrangel).STAsText() as workaround_within_3 + from q + + -> 7669.10402181435 POINT (-179.3 71.27) GEOMETRYCOLLECTION EMPTY GEOMETRYCOLLECTION EMPTY + + PostGIS knows Within for Geography neither, and the intersection trick gives the same result + + */ + + bg::read_wkt("POLYGON((-178.568604 71.564148,-178.017548 71.449692,-177.833313 71.3461,-177.502838 71.277466 ,-177.439453 71.226929,-177.620026 71.116638,-177.9389 71.037491,-178.8186 70.979965,-179.274445 70.907761,-180 70.9972,179.678314 70.895538,179.272766 70.888596,178.791016 70.7964,178.617737 71.035538,178.872192 71.217484,179.530273 71.4383 ,-180 71.535843 ,-179.628601 71.577194,-179.305298 71.551361,-179.03421 71.597748,-178.568604 71.564148))", wrangel); + + bool within = bg::within(Point(-179.3, 71.27), wrangel); + BOOST_CHECK_EQUAL(within, true); + + within = bg::within(Point(-179.9, 70.95), wrangel); + BOOST_CHECK_EQUAL(within, false); + + within = bg::within(Point(179.9, 70.95), wrangel); + BOOST_CHECK_EQUAL(within, false); + + // Test using great circle mapper + // http://www.gcmap.com/mapui?P=5E52N-9E53N-7E50N-5E52N,7E52.5N,8E51.5N,6E51N + + bg::model::polygon triangle; + bg::read_wkt("POLYGON((5 52,9 53,7 50,5 52))", triangle); + BOOST_CHECK_EQUAL(bg::within(Point(7, 52.5), triangle), true); + BOOST_CHECK_EQUAL(bg::within(Point(8.0, 51.5), triangle), false); + BOOST_CHECK_EQUAL(bg::within(Point(6.0, 51.0), triangle), false); +} int test_main( int , char* [] ) @@ -70,8 +114,12 @@ int test_main( int , char* [] ) test_all >(); test_all >(); + test_spherical > >(); + + #if defined(HAVE_TTMATH) test_all >(); + test_spherical > >(); #endif return 0; diff --git a/test/geometry_test_common.hpp b/test/geometry_test_common.hpp index 9fbe8b26d..0af0a8340 100644 --- a/test/geometry_test_common.hpp +++ b/test/geometry_test_common.hpp @@ -94,6 +94,29 @@ template <> struct string_from_type #endif + +struct geographic_policy +{ + template + static inline CoordinateType apply(CoordinateType const& value) + { + return value; + } +}; + +struct mathematical_policy +{ + template + static inline CoordinateType apply(CoordinateType const& value) + { + return 90 - value; + } + +}; + + + + // For all tests: // - do NOT use "using namespace boost::geometry" to make clear what is Boost.Geometry // - use bg:: as short alias diff --git a/test/strategies/Jamfile.v2 b/test/strategies/Jamfile.v2 index 87e1a74e2..08bfeea89 100644 --- a/test/strategies/Jamfile.v2 +++ b/test/strategies/Jamfile.v2 @@ -14,6 +14,7 @@ test-suite boost-geometry-strategies [ run haversine.cpp ] [ run projected_point.cpp ] [ run pythagoras.cpp ] + [ run spherical_side.cpp ] [ run transformer.cpp ] [ run within.cpp ] ; diff --git a/test/strategies/cross_track.cpp b/test/strategies/cross_track.cpp index f53f6255b..debc5aba8 100644 --- a/test/strategies/cross_track.cpp +++ b/test/strategies/cross_track.cpp @@ -28,9 +28,10 @@ #include +// This test is GIS oriented. -template +template void test_distance( typename bg::coordinate_type::type const& lon1, typename bg::coordinate_type::type const& lat1, @@ -60,9 +61,9 @@ void test_distance( Point p1, p2, p3; - bg::assign_values(p1, lon1, lat1); - bg::assign_values(p2, lon2, lat2); - bg::assign_values(p3, lon3, lat3); + bg::assign_values(p1, lon1, LatitudePolicy::apply(lat1)); + bg::assign_values(p2, lon2, LatitudePolicy::apply(lat2)); + bg::assign_values(p3, lon3, LatitudePolicy::apply(lat3)); strategy_type strategy; @@ -83,29 +84,30 @@ void test_distance( } - -template +template void test_all() { typename bg::coordinate_type::type const average_earth_radius = 6372795.0; // distance (Paris <-> Amsterdam/Barcelona), // with coordinates rounded as below ~87 km - // should be is equal - // to distance (Paris <-> Barcelona/Amsterdam) + // is equal to distance (Paris <-> Barcelona/Amsterdam) typename bg::coordinate_type::type const p_to_ab = 86.798321 * 1000.0; - test_distance(2, 48, 4, 52, 2, 41, average_earth_radius, p_to_ab, 0.1); - test_distance(2, 48, 2, 41, 4, 52, average_earth_radius, p_to_ab, 0.1); + test_distance(2, 48, 4, 52, 2, 41, average_earth_radius, p_to_ab, 0.1); + test_distance(2, 48, 2, 41, 4, 52, average_earth_radius, p_to_ab, 0.1); } int test_main(int, char* []) { - test_all > >(); + test_all >, geographic_policy >(); + + // NYI: haversine for mathematical spherical coordinate systems + // test_all >, mathematical_policya >(); #if defined(HAVE_TTMATH) typedef ttmath::Big<1,4> tt; - //test_all > >(); + //test_all >, geographic_policy>(); #endif return 0; diff --git a/test/strategies/haversine.cpp b/test/strategies/haversine.cpp index 3ad545878..7da720652 100644 --- a/test/strategies/haversine.cpp +++ b/test/strategies/haversine.cpp @@ -34,7 +34,7 @@ double const average_earth_radius = 6372795.0; -template +template struct test_distance { typedef bg::strategy::distance::haversine @@ -60,15 +60,15 @@ struct test_distance haversine_type strategy(radius); Point p1, p2; - bg::assign_values(p1, lon1, lat1); - bg::assign_values(p2, lon2, lat2); + 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 +template void test_all() { // earth to unit-sphere -> divide by earth circumference, then it is from 0-1, @@ -77,19 +77,19 @@ void test_all() // ~ Amsterdam/Paris, 467 kilometers double const a_p = 467.2704 * 1000.0; - test_distance::test(4, 52, 2, 48, average_earth_radius, a_p, 1.0); - test_distance::test(2, 48, 4, 52, average_earth_radius, a_p, 1.0); - test_distance::test(4, 52, 2, 48, 1.0, a_p * e2u, 0.001); + test_distance::test(4, 52, 2, 48, average_earth_radius, a_p, 1.0); + test_distance::test(2, 48, 4, 52, average_earth_radius, a_p, 1.0); + test_distance::test(4, 52, 2, 48, 1.0, a_p * e2u, 0.001); // ~ Amsterdam/Barcelona double const a_b = 1232.9065 * 1000.0; - test_distance::test(4, 52, 2, 41, average_earth_radius, a_b, 1.0); - test_distance::test(2, 41, 4, 52, average_earth_radius, a_b, 1.0); - test_distance::test(4, 52, 2, 41, 1.0, a_b * e2u, 0.001); + test_distance::test(4, 52, 2, 41, average_earth_radius, a_b, 1.0); + test_distance::test(2, 41, 4, 52, average_earth_radius, a_b, 1.0); + test_distance::test(4, 52, 2, 41, 1.0, a_b * e2u, 0.001); } -template +template void test_services() { namespace bgsd = bg::strategy::distance; @@ -244,9 +244,12 @@ double time_normal(int n) int test_main(int, char* []) { - test_all > >(); - test_all > >(); - test_all > >(); + test_all >, geographic_policy>(); + test_all >, geographic_policy>(); + test_all >, geographic_policy>(); + + // NYI: haversine for mathematical spherical coordinate systems + // test_all >, mathematical_policy>(); //double t1 = time_sqrt(20000); //double t2 = time_normal(20000); @@ -255,15 +258,16 @@ int test_main(int, char* []) #if defined(HAVE_TTMATH) typedef ttmath::Big<1,4> tt; - test_all > >(); + test_all >, geographic_policy>(); #endif test_services < - bg::model::point >, - bg::model::point >, - double + bg::model::point >, + bg::model::point >, + double, + geographic_policy >(); return 0; diff --git a/test/strategies/segment_intersection.cpp b/test/strategies/segment_intersection.cpp index d68b6ee5c..5dfe3ed0f 100644 --- a/test/strategies/segment_intersection.cpp +++ b/test/strategies/segment_intersection.cpp @@ -91,7 +91,7 @@ static void test_segment_intersection(int caseno, #endif typedef typename bg::coordinate_type

::type coordinate_type; - typedef segment segment_type; + typedef bg::model::referring_segment segment_type; P p1, p2, p3, p4; bg::assign_values(p1, x1, y1); diff --git a/test/strategies/segment_intersection_collinear.cpp b/test/strategies/segment_intersection_collinear.cpp index 3c02f39ad..ea28f9cc5 100644 --- a/test/strategies/segment_intersection_collinear.cpp +++ b/test/strategies/segment_intersection_collinear.cpp @@ -64,7 +64,7 @@ static void test_segment_intersection(std::string const& case_id, //#endif typedef typename bg::coordinate_type

::type coordinate_type; - typedef bg::segment segment_type; + typedef bg::model::referring_segment segment_type; P p1, p2, p3, p4; bg::assign_values(p1, x1, y1); diff --git a/test/strategies/side_by_cross_track.cpp b/test/strategies/side_by_cross_track.cpp deleted file mode 100644 index e02f3267d..000000000 --- a/test/strategies/side_by_cross_track.cpp +++ /dev/null @@ -1,54 +0,0 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) -// Unit Test - -// Copyright (c) 2007-2011 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) - - -#include - -#include - -#include - -#include - -#include -#include - - -template -void test_side(double lon1, double lat1, - double lon2, double lat2, - double lon3, double lat3, - int expected) -{ - typedef bg::strategy::side::side_by_cross_track strategy; - - Point p1, p2, p3; - bg::assign_values(p1, lon1, lat1); - bg::assign_values(p2, lon2, lat2); - bg::assign_values(p3, lon3, lat3); - int s = strategy::apply(p1, p2, p3); - -} - -template -void test_all() -{ - test_side(2.0, 48.0, 4.0, 52.0, 2.0, 41.0, 1); - test_side(2.0, 48.0, 2.0, 41.0, 4.0, 52.0, -1); -} - -int test_main(int, char* []) -{ - test_all > >(); - - double a = 0; - double b = sin(a); - - return 0; -} diff --git a/test/strategies/spherical_side.cpp b/test/strategies/spherical_side.cpp new file mode 100644 index 000000000..5410c6992 --- /dev/null +++ b/test/strategies/spherical_side.cpp @@ -0,0 +1,142 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2007-2011 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) + + +#include + + +#include + + +#include +//#include +#include +#include + +#include + +#include +#include + + +namespace boost { namespace geometry { + +template +static inline Vector create_vector(Point1 const& p1, Point2 const& p2) +{ + Vector v; + convert(p1, v); + subtract_point(v, p2); + return v; +} + +}} + +inline char side_char(int side) +{ + return side == 1 ? 'L' + : side == -1 ? 'R' + : '-' + ; +} + +template +void test_side1(std::string const& case_id, Point const& p1, Point const& p2, Point const& p3, + int expected, int expected_cartesian) +{ + // std::cout << case_id << ": "; + //int s = bg::strategy::side::side_via_plane<>::apply(p1, p2, p3); + int side_ssf = bg::strategy::side::spherical_side_formula<>::apply(p1, p2, p3); + //int side2 = bg::strategy::side::side_via_plane<>::apply(p1, p2, p3); + int side_ct = bg::strategy::side::side_by_cross_track<>::apply(p1, p2, p3); + + typedef bg::strategy::side::services::default_strategy::type cartesian_strategy; + int side_cart = cartesian_strategy::apply(p1, p2, p3); + + + BOOST_CHECK_EQUAL(side_ssf, expected); + BOOST_CHECK_EQUAL(side_ct, expected); + BOOST_CHECK_EQUAL(side_cart, expected_cartesian); + /* + std::cout + << "exp: " << side_char(expected) + << " ssf: " << side_char(side1) + << " pln: " << side_char(side2) + << " ct: " << side_char(side3) + //<< " def: " << side_char(side4) + << " cart: " << side_char(side5) + << std::endl; + */ +} + +template +void test_side(std::string const& case_id, Point const& p1, Point const& p2, Point const& p3, + int expected, int expected_cartesian = -999) +{ + if (expected_cartesian == -999) + { + expected_cartesian = expected; + } + test_side1(case_id, p1, p2, p3, expected, expected_cartesian); + test_side1(case_id, p2, p1, p3, -expected, -expected_cartesian); +} + + +template +void test_all() +{ + typedef std::pair pair; + + Point amsterdam(5.9, 52.4); + Point barcelona(2.0, 41.0); + Point paris(2.0, 48.0); + Point milan(7.0, 45.0); + + //goto wrong; + + test_side("bp-m", barcelona, paris, milan, -1); + test_side("bm-p", barcelona, milan, paris, 1); + test_side("mp-b", milan, paris, barcelona, 1); + + test_side("am-p", amsterdam, milan, paris, -1); + test_side("pm-a", paris, milan, amsterdam, 1); + + // http://www.gcmap.com/mapui?P=30N+10E-50N+50E,39N+30E + Point gcmap_p1(10.0, 30.0); + Point gcmap_p2(50.0, 50.0); + test_side("blog1", gcmap_p1, gcmap_p2, Point(30.0, 41.0), -1, 1); + test_side("blog1", gcmap_p1, gcmap_p2, Point(30.0, 42.0), -1, 1); + test_side("blog1", gcmap_p1, gcmap_p2, Point(30.0, 43.0), -1, 1); + test_side("blog1", gcmap_p1, gcmap_p2, Point(30.0, 44.0), 1); + + // http://www.gcmap.com/mapui?P=50N+80E-60N+50W,65N+30E + Point gcmap_np1(80.0, 50.0); + Point gcmap_np2(-50.0, 60.0); + // http://www.gcmap.com/mapui?P=50N+140E-60N+10E,65N+30E + //Point gcmap_np1(140.0, 50.0); + //Point gcmap_np2(10.0, 60.0); + //test_side(gcmap_np1, gcmap_np2, gcmap_np, 1); + test_side("40", gcmap_np1, gcmap_np2, Point(30.0, 60.0), 1, -1); + test_side("45", gcmap_np1, gcmap_np2, Point(30.0, 65.0), 1, -1); + test_side("70", gcmap_np1, gcmap_np2, Point(30.0, 70.0), 1, -1); + test_side("75", gcmap_np1, gcmap_np2, Point(30.0, 75.0), -1); +} + +int test_main(int, char* []) +{ + test_all > >(); + test_all > >(); + +#if defined(HAVE_TTMATH) + typedef ttmath::Big<1,4> tt; + test_all > >(); +#endif + + return 0; +} diff --git a/test/strategies/side_by_cross_track.vcproj b/test/strategies/spherical_side.vcproj similarity index 92% rename from test/strategies/side_by_cross_track.vcproj rename to test/strategies/spherical_side.vcproj index e5f1a9359..70dd8689c 100644 --- a/test/strategies/side_by_cross_track.vcproj +++ b/test/strategies/spherical_side.vcproj @@ -2,9 +2,9 @@ @@ -18,7 +18,7 @@ diff --git a/test/strategies/strategies_tests.sln b/test/strategies/strategies_tests.sln index 6120639f4..226baff48 100644 --- a/test/strategies/strategies_tests.sln +++ b/test/strategies/strategies_tests.sln @@ -14,7 +14,7 @@ Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "projected_point", "projecte EndProject Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "segment_intersection_collinear", "segment_intersection_collinear.vcproj", "{2D0CB6D3-6ABC-4119-A235-66E6065A279E}" EndProject -Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "side_by_cross_track", "side_by_cross_track.vcproj", "{ADBE38D8-1828-48A2-BBA1-81F50B53C67C}" +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "spherical_side", "spherical_side.vcproj", "{ADBE38D8-1828-48A2-BBA1-81F50B53C67C}" EndProject Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "within", "within.vcproj", "{AB13D2AC-FD34-4DE4-BD8E-4D463050E5DD}" EndProject