From 5cdae7c693b6c5290bfcb9db79d02af0cecfa97b Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 2 May 2017 18:44:26 +0300 Subject: [PATCH] [disjoint] [vertex_longitude] Reviewing cases of disjoint, create new tests, fix corner cases --- .../detail/disjoint/segment_box.hpp | 155 ++++---- .../geometry/formulas/vertex_longitude.hpp | 48 ++- .../disjoint/disjoint_seg_box.cpp | 130 ++++++- test/formulas/vertex_longitude.cpp | 24 +- test/formulas/vertex_longitude_cases.hpp | 342 +++++++++++------- 5 files changed, 463 insertions(+), 236 deletions(-) diff --git a/include/boost/geometry/algorithms/detail/disjoint/segment_box.hpp b/include/boost/geometry/algorithms/detail/disjoint/segment_box.hpp index b2ad15f3e..04e1fae3c 100644 --- a/include/boost/geometry/algorithms/detail/disjoint/segment_box.hpp +++ b/include/boost/geometry/algorithms/detail/disjoint/segment_box.hpp @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -36,6 +37,7 @@ #include #include +#include //TODO:remove if not needed namespace boost { namespace geometry { @@ -77,14 +79,15 @@ public: geometry::detail::assign_point_from_index<0>(segment, p0); geometry::detail::assign_point_from_index<1>(segment, p1); - // Test simple case of intersection first + // Simplest cases first + // Case 1: if box contains one of segment's endpoints then they are not disjoint if (! disjoint_point_box(p0, box) || ! disjoint_point_box(p1, box)) { return false; } - // Test intersection by comparing angles + // Case 2: disjoint if bounding boxes are disjoint typedef typename coordinate_type::type CT; @@ -103,34 +106,9 @@ public: swap(lon1, lat1, lon2, lat2); } - CT alp1, a_b0, a_b1, a_b2, a_b3; - - CT b_lon_min = geometry::get_as_radian(box); - CT b_lat_min = geometry::get_as_radian(box); - CT b_lon_max = geometry::get_as_radian(box); - CT b_lat_max = geometry::get_as_radian(box); + CT alp1; azimuth_strategy.apply(lon1, lat1, lon2, lat2, alp1); - azimuth_strategy.apply(lon1, lat1, b_lon_min, b_lat_min, a_b0); - azimuth_strategy.apply(lon1, lat1, b_lon_max, b_lat_min, a_b1); - azimuth_strategy.apply(lon1, lat1, b_lon_min, b_lat_max, a_b2); - azimuth_strategy.apply(lon1, lat1, b_lon_max, b_lat_max, a_b3); - - bool b0 = alp1 > a_b0; - bool b1 = alp1 > a_b1; - bool b2 = alp1 > a_b2; - bool b3 = alp1 > a_b3; - - // if the box is above (below) the segment in northern (southern) - // hemisphere respectively then there is not intersection - if (!(b0 && b1 && b2 && b3) && (b0 || b1 || b2 || b3)) - { - return false; - } - - // The only case not covered by the angle test above - // Test intersection by testing if the vertex of the geodesic segment - // is covered by the box geometry::model::box box_seg; @@ -141,58 +119,93 @@ public: box_seg, azimuth_strategy, alp1); + if (disjoint_box_box(box, box_seg)) + { + return true; + } + + // Case 3: test intersection by comparing angles + + CT a_b0, a_b1, a_b2, a_b3; + + CT b_lon_min = geometry::get_as_radian(box); + CT b_lat_min = geometry::get_as_radian(box); + CT b_lon_max = geometry::get_as_radian(box); + CT b_lat_max = geometry::get_as_radian(box); + + azimuth_strategy.apply(lon1, lat1, b_lon_min, b_lat_min, a_b0); + azimuth_strategy.apply(lon1, lat1, b_lon_max, b_lat_min, a_b1); + azimuth_strategy.apply(lon1, lat1, b_lon_min, b_lat_max, a_b2); + azimuth_strategy.apply(lon1, lat1, b_lon_max, b_lat_max, a_b3); + + bool b0 = alp1 > a_b0; + bool b1 = alp1 > a_b1; + bool b2 = alp1 > a_b2; + bool b3 = alp1 > a_b3; + + // if not all box points on the same side of the segment then + // there is an intersection + if (!(b0 && b1 && b2 && b3) && (b0 || b1 || b2 || b3)) + { + return false; + } + + // Case 4: The only intersection case not covered above is when all four + // points of the box are above (below) the segment in northern (southern) + // hemisphere. Then we have to compute the vertex of the segment CT vertex_lat; - if (lat1 < CT(0) && lat2 < CT(0)) + CT lat_sum = lat1 + lat2; + + if ((b0 && b1 && b2 && b3 && lat_sum > CT(0)) + || (!(b0 && b1 && b2 && b3) && lat_sum < CT(0))) { - vertex_lat = geometry::get_as_radian(box); - } - if (lat1 > CT(0) && lat2 > CT(0)) - { - vertex_lat = geometry::get_as_radian(box); - } - if (lat1 > CT(0) && lat2 < CT(0)) - { - if (lat1 > -lat2) + CT b_lat_below; //latitude of box closest to equator + + if (lat_sum > CT(0)) { - vertex_lat = geometry::get_as_radian(box); + vertex_lat = geometry::get_as_radian(box_seg); + b_lat_below = b_lat_min; } else { - vertex_lat = geometry::get_as_radian(box); + vertex_lat = geometry::get_as_radian(box_seg); + b_lat_below = b_lat_max; } - } - if (lat1 < CT(0) && lat2 > CT(0)) - { - if (-lat1 < lat2) + + //TODO: computing the spherical longitude should suffice for this test (?) + CT vertex_lon = geometry::formula::vertex_longitude::apply(lon1, lat1, + lon2, lat2, + vertex_lat, + alp1, + azimuth_strategy); + + segment_point_type p_vertex_rad; + geometry::set_from_radian<0>(p_vertex_rad, vertex_lon); + geometry::set_from_radian<1>(p_vertex_rad, vertex_lat); + + std::cout << "vertex=" << vertex_lon * math::r2d() + << " , " << vertex_lat * math::r2d() + << std::endl << std::endl; + + + Box box_rad; + geometry::set_from_radian(box_rad, b_lon_min); + geometry::set_from_radian(box_rad, b_lat_min); + geometry::set_from_radian(box_rad, b_lon_max); + geometry::set_from_radian(box_rad, b_lat_max); + + + // Check if the vertex point is within the band defined by the + // minimum and maximum longitude of the box; if yes, then return + // false if the point is above the min latitude of the box; return + // true in all other cases + if (vertex_lon >= b_lon_min && vertex_lon <= b_lon_max + && std::abs(vertex_lat) > std::abs(b_lat_below)) { - vertex_lat = geometry::get_as_radian(box); - } else { - vertex_lat = geometry::get_as_radian(box); + return false; } } - - CT vertex_lon = geometry::formula::vertex_longitude::apply(lon1, lat1, - lon2, lat2, - vertex_lat, - alp1, - azimuth_strategy); - - segment_point_type p_vertex_rad; - geometry::set_from_radian<0>(p_vertex_rad, vertex_lon); - geometry::set_from_radian<1>(p_vertex_rad, vertex_lat); - - //std::cout << "vertex=" << vertex_lon * math::r2d() - // << " , " << vertex_lat * math::r2d() - // << std::endl; - - - Box box_rad; - geometry::set_from_radian(box_rad, b_lon_min); - geometry::set_from_radian(box_rad, b_lat_min); - geometry::set_from_radian(box_rad, b_lon_max); - geometry::set_from_radian(box_rad, b_lat_max); - - return disjoint_point_box(p_vertex_rad, box_rad); + return true; } }; @@ -218,7 +231,7 @@ namespace dispatch template struct disjoint - : detail::disjoint::disjoint_segment_box + : detail::disjoint::disjoint_segment_box {}; diff --git a/include/boost/geometry/formulas/vertex_longitude.hpp b/include/boost/geometry/formulas/vertex_longitude.hpp index 123147083..3c1ae6bea 100644 --- a/include/boost/geometry/formulas/vertex_longitude.hpp +++ b/include/boost/geometry/formulas/vertex_longitude.hpp @@ -220,23 +220,30 @@ public: omg12 += omg1 - omg2; - CT omg13 = vertex_longitude_on_sphere - ::apply(bet1, bet2, bet3, sin(omg12), cos(omg12)); + CT sin_omg12 = sin(omg12); + CT cos_omg12 = cos(omg12); - if (sin(omg12) > 0 && cos(omg12) < 0) + CT omg13 = geometry::formula::vertex_longitude_on_sphere + ::apply(bet1, bet2, bet3, sin_omg12, cos_omg12); + + + if (lat1 * lat2 < CT(0))//different hemispheres { - omg13 = pi - omg13; + if ((lat2-lat1)*lat3 > 0)// ascending segment + { + omg13 = pi - omg13; + } } - std::cout << "bet=" << bet1 * geometry::math::r2d() - << " " << bet2 * geometry::math::r2d() - << " " << bet3 * geometry::math::r2d() - << "\n omg12= " << (omg1 - omg2) * geometry::math::r2d() - << " " << omg12 * geometry::math::r2d() - << "\n omg1= " << omg1 * geometry::math::r2d() - << "\n omg2= " << omg2 * geometry::math::r2d() + std::cout << "bet=" << bet1 * geometry::math::r2d() + << " " << bet2 * geometry::math::r2d() + << " " << bet3 * geometry::math::r2d() + << "\n omg12= " << (omg1 - omg2) * geometry::math::r2d() + << " " << omg12 * geometry::math::r2d() + << "\n omg1= " << omg1 * geometry::math::r2d() + << "\n omg2= " << omg2 * geometry::math::r2d() << "\n omg13(rad)= " << omg13 - << "\n omg13= " << omg13 * geometry::math::r2d() + << "\n omg13= " << omg13 * geometry::math::r2d() << "\n"; // Second, compute the ellipsoidal longitude @@ -278,7 +285,7 @@ public: + C32 * (sin4_sig3 - sin4_sig1)); - std::cout << "correction=" << f * sin_alp0 * I3 * geometry::math::r2d() << std::endl; + std::cout << "correction=" << f * sin_alp0 * I3 * geometry::math::r2d() << std::endl; std::cout << "I3=" << I3 << std::endl; int sign = CT(1); @@ -349,6 +356,7 @@ struct compute_vertex_lon }; // Vertex longitude interface +// Assume that lon1 < lon2 and lon1, lon2 \in (-pi,pi] template class vertex_longitude @@ -423,14 +431,24 @@ public : azimuth_strategy); CT pi = math::pi(); - CT vertex_lon = std::fmod(lon1 + dlon + pi, 2 * pi) - pi; + CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi);// - pi;; + if(vertex_lat < CT(0)) { - vertex_lon += geometry::math::pi(); + vertex_lon -= geometry::math::pi(); } + if (std::abs(lon1 - lon2) > pi) + { + vertex_lon -= pi; + } + + vertex_lon = std::fmod(vertex_lon, 2 * pi); + + //geometry::math::normalize_longitude(vertex_lon); + return vertex_lon; diff --git a/test/algorithms/relational_operations/disjoint/disjoint_seg_box.cpp b/test/algorithms/relational_operations/disjoint/disjoint_seg_box.cpp index b2f79340c..06f67d5ff 100644 --- a/test/algorithms/relational_operations/disjoint/disjoint_seg_box.cpp +++ b/test/algorithms/relational_operations/disjoint/disjoint_seg_box.cpp @@ -11,6 +11,8 @@ #define BOOST_GEOMETRY_TEST_DEBUG +#include + #include #include #include @@ -30,6 +32,7 @@ #include "test_disjoint_seg_box.hpp" + namespace bg = boost::geometry; //Tests for disjoint(point, box), disjoint(box, box) and disjoint(segment, box) @@ -73,7 +76,6 @@ void disjoint_tests_1() test_disjoint, bg::model::segment

>("BOX(1 1,3 3)", "SEGMENT(0 0, 4 0)", true); - test_disjoint, bg::model::segment

>("BOX(1 1,3 3)", "SEGMENT(2 2, 4 4)", false); @@ -122,25 +124,25 @@ template void disjoint_tests_with_strategy(bool expected_result) { bg::strategy::disjoint::segment_box_geographic - < - bg::strategy::andoyer, - bg::srs::spheroid, - CT - > geographic_andoyer; + < + bg::strategy::andoyer, + bg::srs::spheroid, + CT + > geographic_andoyer; bg::strategy::disjoint::segment_box_geographic - < - bg::strategy::thomas, - bg::srs::spheroid, - CT - > geographic_thomas; + < + bg::strategy::thomas, + bg::srs::spheroid, + CT + > geographic_thomas; bg::strategy::disjoint::segment_box_geographic - < - bg::strategy::vincenty, - bg::srs::spheroid, - CT - > geographic_vincenty; + < + bg::strategy::vincenty, + bg::srs::spheroid, + CT + > geographic_vincenty; test_disjoint_strategy, bg::model::segment

> ("BOX(1 1,3 3)", "SEGMENT(1 0.999, 10 0.999)", @@ -153,6 +155,99 @@ void disjoint_tests_with_strategy(bool expected_result) expected_result, geographic_vincenty); } +template +void disjoint_tests_sph_geo() +{ + //Case A: box intersects without containing the vertex + test_disjoint, bg::model::segment

>("BOX(0 6, 120 7)", + "SEGMENT(0 5, 120 5)", + false); + test_disjoint, bg::model::segment

>("BOX(0 -6, 120 -7)", + "SEGMENT(0 -5, 120 -5)", + false); + + //Case B: box intersects and contains the vertex + test_disjoint, bg::model::segment

>("BOX(0 9, 120 10)", + "SEGMENT(0 5, 120 5)", + false); + test_disjoint, bg::model::segment

>("BOX(0 -10, 120 -9)", + "SEGMENT(0 -5, 120 -5)", + false); + + //Case C: bounding boxes disjoint + test_disjoint, bg::model::segment

>("BOX(0 10, 10 20)", + "SEGMENT(0 5, 120 5)", + true); + test_disjoint, bg::model::segment

>("BOX(0 -20, 10 -10)", + "SEGMENT(0 -5, 120 -5)", + true); + + //Case D: bounding boxes intersect but box segment are disjoint + test_disjoint, bg::model::segment

>("BOX(0 9, 0.1 20)", + "SEGMENT(0 5, 120 5)", + true); + test_disjoint, bg::model::segment

>("BOX(0 -20, 0.1 -9)", + "SEGMENT(0 -5, 120 -5)", + true); + + //Case E: geodesic intersects box but box segment are disjoint + test_disjoint, bg::model::segment

>("BOX(121 0, 122 10)", + "SEGMENT(0 5, 120 5)", + true); + test_disjoint, bg::model::segment

>("BOX(121 -10, 122 0)", + "SEGMENT(0 -5, 120 -5)", + true); + + //Case F: segment crosses box + test_disjoint, bg::model::segment

>("BOX(100 0, 110 20)", + "SEGMENT(0 5, 120 5)", + false); + test_disjoint, bg::model::segment

>("BOX(100 -20, 110 0)", + "SEGMENT(0 -5, 120 -5)", + false); + + //Case G: box contains one segment endpoint + test_disjoint, bg::model::segment

>("BOX(110 0, 130 10)", + "SEGMENT(0 5, 120 5)", + false); + test_disjoint, bg::model::segment

>("BOX(110 -10, 130 0)", + "SEGMENT(0 -5, 120 -5)", + false); + + //Case H: box below segment + test_disjoint, bg::model::segment

>("BOX(50 0, 70 6)", + "SEGMENT(0 5, 120 5)", + true); + test_disjoint, bg::model::segment

>("BOX(50 -6, 70 0)", + "SEGMENT(0 -5, 120 -5)", + true); + + //Case I: box contains segment + test_disjoint, bg::model::segment

>("BOX(-10 0, 130 10)", + "SEGMENT(0 5, 120 5)", + false); + test_disjoint, bg::model::segment

>("BOX(-10 -10, 130 0)", + "SEGMENT(0 -5, 120 -5)", + false); + + //ascending segment + test_disjoint, bg::model::segment

>("BOX(0 10, 120 10.1)", + "SEGMENT(0 5, 120 5.1)", + false); + //descending segment + test_disjoint, bg::model::segment

>("BOX(0 9.8, 120 10)", + "SEGMENT(0 5, 120 4.9)", + false); + //ascending segment both hemispheres + test_disjoint, bg::model::segment

>("BOX(100 5, 120 6)", + "SEGMENT(0 -1, 120 4.9)", + false); + //descending segment both hemispheres + test_disjoint, bg::model::segment

>("BOX(0 5, 20 6)", + "SEGMENT(0 4.9, 120 -1)", + false); +} + template void test_all() { @@ -179,6 +274,9 @@ void test_all() disjoint_tests_4(false); disjoint_tests_with_strategy(false); + + disjoint_tests_sph_geo(); + disjoint_tests_sph_geo(); } diff --git a/test/formulas/vertex_longitude.cpp b/test/formulas/vertex_longitude.cpp index 9fd1484f9..bbd945355 100644 --- a/test/formulas/vertex_longitude.cpp +++ b/test/formulas/vertex_longitude.cpp @@ -13,7 +13,7 @@ #include "test_formula.hpp" #include "vertex_longitude_cases.hpp" -#include +#include //TODO:remove if not needed #include #include @@ -25,6 +25,8 @@ #include +#define BOOST_GEOMETRY_TEST_DEBUG + namespace bg = boost::geometry; template @@ -66,7 +68,6 @@ CT test_vrt_lon_sph(CT lon1r, vertex_lat, a1, azimuth); - } template @@ -125,12 +126,22 @@ CT test_vrt_lon_geo(CT lon1r, void test_all(expected_results const& results) { double const d2r = bg::math::d2r(); + double const pi = bg::math::pi(); double lon1r = results.p1.lon * d2r; double lat1r = results.p1.lat * d2r; double lon2r = results.p2.lon * d2r; double lat2r = results.p2.lat * d2r; - +/* + if (lon1r > pi) + { + lon1r -= 2 * pi; + } + if (lon2r > pi) + { + lon2r -= 2 * pi; + } +*/ if(lon1r > lon2r) { std::swap(lon1r, lon2r); @@ -151,8 +162,13 @@ void test_all(expected_results const& results) std::cout << res_vi * bg::math::r2d() << "," << std::endl; std::cout << res_sh * bg::math::r2d() << std::endl<< std::endl; + bg::math::normalize_longitude(res_an); + bg::math::normalize_longitude(res_th); + bg::math::normalize_longitude(res_vi); + bg::math::normalize_longitude(res_sh); + check_one(res_an, results.andoyer * d2r, res_vi, 0.001); - check_one(res_th, results.thomas * d2r, res_vi, 0.00001); + check_one(res_th, results.thomas * d2r, res_vi, 0.01);//in some tests thomas gives low accuracy check_one(res_vi, results.vincenty * d2r, res_vi, 0.0000001); check_one(res_sh, results.spherical * d2r, res_vi, 1); } diff --git a/test/formulas/vertex_longitude_cases.hpp b/test/formulas/vertex_longitude_cases.hpp index a116375a1..d42d77c38 100644 --- a/test/formulas/vertex_longitude_cases.hpp +++ b/test/formulas/vertex_longitude_cases.hpp @@ -35,8 +35,7 @@ struct expected_results }; expected_results expected[] = -{ - +{ { //ascenting segments (wrt pole) { 1, 1 },{ 100, 2 }, 66.25553538, @@ -121,6 +120,30 @@ expected_results expected[] = 1.935976226, 1.93597805, 1.779116148 + },{ + { 3, 5 },{ 150, 0}, + 60.29182988, + 60.29785309, + 60.29785255, + 60 + },{ + { 3, 5 },{ 150, 0.5}, + 63.11344576, + 63.11900045, + 63.11899891, + 62.87000766 + },{ + { 3, 5 },{ 150, 1}, + 65.51880171, + 65.52391866, + 65.52391623, + 65.31813729 + },{ + { 3, 5 },{ 150, 5}, + 76.49727275, + 76.50000657, + 76.5, + 76.5 },{ //segments parallel to equator { 0, 1 },{ 4, 1 }, 1.999999973, @@ -140,29 +163,35 @@ expected_results expected[] = 30, 30 },{ - { 0, 1 },{ 180, 1 }, - 90, - 90, - 90, - 90 - },/*{ { 0, 1 },{ 90, 1 }, 44.99960266, - 45.22272756, + 45.22272756,//thomas low accuracy 45, 45 - },{ - { 0, 1 },{ 270, 1 }, - 90, - 90, - 90, - 90 - },*/{ + },{ { 0, 1 },{ 120, 1 }, 59.99878311, 59.99999778, 60, 60 + },{ + { 0, 1 },{ 180, 1 }, + 90, + 90, + 90, + 90 + },{ + { 0, 1 },{ 270, 1 }, + -44.99960266, + -45.08931472,//thomas low accuracy + -45, + -45 + },{ + { 0, 1 },{ 290, 1 }, + -34.9998314, + -34.99999868, + -35, + -35 },{ { 0, 1 },{ 150, 1 }, 74.99598515, @@ -235,7 +264,38 @@ expected_results expected[] = 24.99999901, 24.99999995, 25 - },{ //Both hemispheres, vertex on the northern + },{//Both hemispheres, vertex on the northern + //A desc vertex north + { 3, 5 },{ 150, -3}, + 27.357069, + 27.36422922, + 27.36423549, + 26.74999989 + },{//B asc vertex north + { 3, -3 },{ 150, 5}, + 125.6403436, + 125.6357677, + 125.6357659, + 126.2500001 + },{//C desc vertex south + { 3, -5 },{ 150, 3}, + 27.3570679, + 27.36422812, + 27.36423439, + 26.74999989 + },{//D asc vertex south + { 3, 3 },{ 150, -5}, + 125.6403423, + 125.6357664, + 125.6357645, + 126.2500001 + },{//E asc vertex south + { 3, 3 },{ 184, -5}, + -88.00743796, + -88.0660268, + -88.0558747, + -88.49315894 + },{ { 3, 5 },{ 150, -3.5}, 17.96722293, 17.97322407, @@ -247,6 +307,30 @@ expected_results expected[] = 52.97759463, 52.9775964, 52.56504545 + },{ //Both hemispheres, vertex on the southern + { 3, 3},{ 5, -5}, + 5, + 5, + 5, + 5 + },{ + { 3, -5 },{ 150, 1}, //symmetric to { 3, 5 },{ 150, -1} + 52.97060093, + 52.97759176, + 52.97759353, + 52.56504545 + },{// fix p1 lon, lat and p2 lon and vary p2 lat + { 3, 5 },{ 150, 1}, + 65.51880171, + 65.52391866, + 65.52391623, + 65.31813729 + },{ + { 3, 5 },{ 150, 0}, + 60.29182988, + 60.29785309, + 60.29785255, + 60 },{ { 3, 5 },{ 150, -0.1}, 59.66911673, @@ -254,121 +338,119 @@ expected_results expected[] = 59.67523616, 59.36690727 },{ - { 3, 5 },{ 150, 0}, - 60.29182988, - 60.29785309, - 60.29785255, - 60 + { 3, 5 },{ 150, -1}, + 52.9706038, + 52.97759463, + 52.9775964, + 52.56504545 },{ - { 3, 5 },{ 150, 1}, - 65.51880171, - 65.52391866, - 65.52391623, - 65.31813729 - },{ //Both hemispheres, vertex on the southern - { 3, 3},{ 5, -5}, - 5, - 5, - 5, - 5 - },{//symmetric to { 3, 5 },{ 150, -1} - { 3, -5 },{ 150, 1}, - 52.97060093, - 52.97759176, - 52.97759353, - 52.56504545 - },{//A desc north - { 3, 5 },{ 150, -3}, - 27.357069, - 27.36422922, - 27.36423549, - 26.74999989 - },{//B asc north - { 3, -3 },{ 150, 5}, - 125.6403436, - 125.6357677, - 125.6357659, - 126.2500001 - },{//C desc south - { 150, 3},{ 3, -5 }, - 27.3570679, - 27.36422812, - 27.36423439, - 26.74999989 - },{//D asc south - { 3, 3 },{ 150, -5}, - 125.6403423, - 125.6357664, - 125.6357645, - 126.2500001 + { 3, 5 },{ 150, -4.15}, + 4.481947557, + 4.485467841, + 4.485473295, + 3.981178967 },{ - { 0, 0.1 },{ 120, 0.05}, - 49.22767004, - 49.22917982, - 49.22918305, - 49.10659455 - },{ - { 0, 0.1 },{ 120, 0}, - 30.30011577, - 30.30174782, - 30.30175235, - 30 - },{ - { 0, 0.1 },{ 120, -0.5}, - 100.5261824, - 100.5257736, - 100.5257733, - 100.8931067 - },{ - { 3, 5 },{ 150, 0}, - 60.29182988, - 60.29785309, - 60.29785255, - 60 - },{ - { 3, 5 },{ 150, 0.5}, - 63.11344576, - 63.11900045, - 63.11899891, - 62.87000766 - },{ - { 3, 5 },{ 150, 1}, - 65.51880171, - 65.52391866, - 65.52391623, - 65.31813729 - },{ - { 3, 5 },{ 150, 5}, - 76.49727275, - 76.50000657, - 76.5, - 76.5 - },{//symmetry of geodesics (1) - { 0, 5 },{ 30, 5.5}, - 25.06431998, - 25.0644277, - 25.06442787, - 25.13253724 - },{//symmetry of geodesics (2) - { 0, 5.5 },{ 30, 5}, - 4.935667094, - 4.935571216, - 4.93557213, - 4.867462762 - },{//symmetry of geodesics (3) - { 0, -5 },{ 30, -5.5}, - 25.06431885, - 25.06442657, - 25.06442674, - 25.13253724 - },{//symmetry of geodesics (4) - { 0, -5.5 },{ 30, -5}, - 4.935666841, - 4.935570963, - 4.935571877, - 4.867462762 - } - + { 3, 5 },{ 150, -4.2}, + 3, + 3, + 3, + 3 + },{//symmetry of geodesics: + // (i) case A same as C and B same as D + // (ii) longitude diff between vertex and p2 in A, C equals + // longitude diff between vertex and p1 in B, D by symmetry + // case (A) + { 0, 5 },{ 30, 5.5}, + 25.06431998, + 25.0644277, + 25.06442787, + 25.13253724 + },{// case (B) + { 0, 5.5 },{ 30, 5}, + 4.935667094, + 4.935571216, + 4.93557213, + 4.867462762 + },{// case (C) + { 0, -5 },{ 30, -5.5}, + 25.06431885, + 25.06442657, + 25.06442674, + 25.13253724 + },{// case (D) + { 0, -5.5 },{ 30, -5}, + 4.935666841, + 4.935570963, + 4.935571877, + 4.867462762 + },{//crossing meridian + { -10, 1 },{ 50, 1.1}, + 24.68113946, + 24.68127641, + 24.68127733, + 24.71605263 + },{ + { 350, 1 },{ 50, 1.1}, + 24.68113946, + 24.68127641, + 24.68127733, + 24.71605263 + },{//crossing antimeridian + { 130, 1 },{ 190, 1.1}, + 164.6811395, + 164.6812764, + 164.6812773, + 164.7160526 + },{ + { 130, 1 },{ -170, 1.1}, + 164.6811395, + 164.6812764, + 164.6812773, + 164.7160526 + },{//crossing meridian both hemispheres + { -10, -5 },{ 150, 1}, + 55.61285835, + 55.62727853, + 55.62725182, + 55.19943725 + },{ + { 350, -5 },{ 150, 1}, + 55.6243632, + 55.6272619, + 55.627257, + 55.1994373 + },{//crossing anti-meridian both hemispheres + { 90, -5 },{ 210, 1}, + 109.4997596, + 109.5011987, + 109.5012031, + 109.1354089 + },{ + { 90, -5 },{ -150, 1}, + 109.4997596, + 109.5011987, + 109.5012031, + 109.1354089 + },{ + { -150, -5 },{ 90, 1}, + -169.4997596, + -169.5011987, + -169.5012031, + -169.1354089 + },{ + { 90, 1 },{ 210, -5}, + -169.5008004, + -169.5012037, + -169.501204, + -169.1354089 + },{ + { 0, 1 },{ 120, -5}, + 100.4991996, + 100.4987963, + 100.498796, + 100.8645911 + } + /**/ }; size_t const expected_size = sizeof(expected) / sizeof(expected_results);