1 // Boost.Geometry
2 // Unit Test
3
4 // Copyright (c) 2016-2019 Oracle and/or its affiliates.
5
6 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
8
9 // Use, modification and distribution is subject to the Boost Software License,
10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
11 // http://www.boost.org/LICENSE_1_0.txt)
12
13 #include "test_formula.hpp"
14 #include "intersection_cases.hpp"
15
16 #include <boost/geometry/formulas/andoyer_inverse.hpp>
17 #include <boost/geometry/formulas/geographic.hpp>
18 #include <boost/geometry/formulas/gnomonic_intersection.hpp>
19 #include <boost/geometry/formulas/sjoberg_intersection.hpp>
20 #include <boost/geometry/formulas/thomas_direct.hpp>
21 #include <boost/geometry/formulas/thomas_inverse.hpp>
22 #include <boost/geometry/formulas/vincenty_direct.hpp>
23 #include <boost/geometry/formulas/vincenty_inverse.hpp>
24
25 #include <boost/geometry/strategies/geographic/parameters.hpp>
26
27 #include <boost/geometry/srs/spheroid.hpp>
28
check_result(expected_result const & result,expected_result const & expected,expected_result const & reference,double reference_error,bool check_reference_only)29 void check_result(expected_result const& result, expected_result const& expected,
30 expected_result const& reference, double reference_error,
31 bool check_reference_only)
32 {
33 //BOOST_CHECK_MESSAGE((false), "(" << result.lon << " " << result.lat << ") vs (" << expected.lon << " " << expected.lat << ")");
34 check_one(result.lon, expected.lon, reference.lon, reference_error, false, check_reference_only);
35 check_one(result.lat, expected.lat, reference.lat, reference_error, false, check_reference_only);
36 }
37
test_formulas(expected_results const & results,bool check_reference_only)38 void test_formulas(expected_results const& results, bool check_reference_only)
39 {
40 // reference result
41 if (results.sjoberg_vincenty.lon == ND)
42 {
43 return;
44 }
45
46 double const d2r = bg::math::d2r<double>();
47 double const r2d = bg::math::r2d<double>();
48
49 double lona1r = results.p1.lon * d2r;
50 double lata1r = results.p1.lat * d2r;
51 double lona2r = results.p2.lon * d2r;
52 double lata2r = results.p2.lat * d2r;
53 double lonb1r = results.q1.lon * d2r;
54 double latb1r = results.q1.lat * d2r;
55 double lonb2r = results.q2.lon * d2r;
56 double latb2r = results.q2.lat * d2r;
57
58 expected_result result;
59
60 // WGS84
61 bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
62
63 if (results.gnomonic_vincenty.lon != ND)
64 {
65 bg::formula::gnomonic_intersection<double, bg::formula::vincenty_inverse, bg::formula::vincenty_direct>
66 ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
67 result.lon *= r2d;
68 result.lat *= r2d;
69 check_result(result, results.gnomonic_vincenty, results.sjoberg_vincenty, 0.00000001, check_reference_only);
70 }
71
72 if (results.gnomonic_thomas.lon != ND)
73 {
74 bg::formula::gnomonic_intersection<double, bg::strategy::thomas::inverse, bg::strategy::thomas::direct>
75 ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
76 result.lon *= r2d;
77 result.lat *= r2d;
78 check_result(result, results.gnomonic_thomas, results.sjoberg_vincenty, 0.0000001, check_reference_only);
79 }
80
81 if (results.sjoberg_vincenty.lon != ND)
82 {
83 bg::formula::sjoberg_intersection<double, bg::formula::vincenty_inverse, 4>
84 ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
85 result.lon *= r2d;
86 result.lat *= r2d;
87 check_result(result, results.sjoberg_vincenty, results.sjoberg_vincenty, 0.00000001, check_reference_only);
88 }
89
90 if (results.sjoberg_thomas.lon != ND)
91 {
92 bg::formula::sjoberg_intersection<double, bg::formula::thomas_inverse, 2>
93 ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
94 result.lon *= r2d;
95 result.lat *= r2d;
96 check_result(result, results.sjoberg_thomas, results.sjoberg_vincenty, 0.0000001, check_reference_only);
97 }
98
99 if (results.sjoberg_andoyer.lon != ND)
100 {
101 bg::formula::sjoberg_intersection<double, bg::formula::andoyer_inverse, 1>
102 ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
103 result.lon *= r2d;
104 result.lat *= r2d;
105 check_result(result, results.sjoberg_andoyer, results.sjoberg_vincenty, 0.0001, check_reference_only);
106 }
107
108 if (results.great_elliptic.lon != ND)
109 {
110 typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree> > point_geo;
111 typedef bg::model::point<double, 3, bg::cs::cartesian> point_3d;
112 point_geo a1(results.p1.lon, results.p1.lat);
113 point_geo a2(results.p2.lon, results.p2.lat);
114 point_geo b1(results.q1.lon, results.q1.lat);
115 point_geo b2(results.q2.lon, results.q2.lat);
116 point_3d a1v = bg::formula::geo_to_cart3d<point_3d>(a1, spheroid);
117 point_3d a2v = bg::formula::geo_to_cart3d<point_3d>(a2, spheroid);
118 point_3d b1v = bg::formula::geo_to_cart3d<point_3d>(b1, spheroid);
119 point_3d b2v = bg::formula::geo_to_cart3d<point_3d>(b2, spheroid);
120 point_3d resv(0, 0, 0);
121 point_geo res(0, 0);
122 bg::formula::great_elliptic_intersection(a1v, a2v, b1v, b2v, resv, spheroid);
123 res = bg::formula::cart3d_to_geo<point_geo>(resv, spheroid);
124 result.lon = bg::get<0>(res);
125 result.lat = bg::get<1>(res);
126 check_result(result, results.great_elliptic, results.sjoberg_vincenty, 0.01, check_reference_only);
127 }
128 }
129
test_4_input_combinations(expected_results const & results,bool check_reference_only)130 void test_4_input_combinations(expected_results const& results, bool check_reference_only)
131 {
132 test_formulas(results, check_reference_only);
133
134 #ifdef BOOST_GEOMETRY_TEST_GEO_INTERSECTION_TEST_SIMILAR
135 {
136 expected_results results_alt = results;
137 std::swap(results_alt.p1, results_alt.p2);
138 test_formulas(results_alt, true);
139 }
140 {
141 expected_results results_alt = results;
142 std::swap(results_alt.q1, results_alt.q2);
143 test_formulas(results_alt, true);
144 }
145 {
146 expected_results results_alt = results;
147 std::swap(results_alt.p1, results_alt.p2);
148 std::swap(results_alt.q1, results_alt.q2);
149 test_formulas(results_alt, true);
150 }
151 #endif
152 }
153
test_all(expected_results const & results)154 void test_all(expected_results const& results)
155 {
156 test_4_input_combinations(results, false);
157
158 #ifdef BOOST_GEOMETRY_TEST_GEO_INTERSECTION_TEST_SIMILAR
159 expected_results results_alt = results;
160 results_alt.p1.lat *= -1;
161 results_alt.p2.lat *= -1;
162 results_alt.q1.lat *= -1;
163 results_alt.q2.lat *= -1;
164 results_alt.gnomonic_vincenty.lat *= -1;
165 results_alt.gnomonic_thomas.lat *= -1;
166 results_alt.sjoberg_vincenty.lat *= -1;
167 results_alt.sjoberg_thomas.lat *= -1;
168 results_alt.sjoberg_andoyer.lat *= -1;
169 results_alt.great_elliptic.lat *= -1;
170 test_4_input_combinations(results_alt, true);
171 #endif
172 }
173
test_bugs()174 void test_bugs()
175 {
176 // https://github.com/boostorg/geometry/issues/612
177 {
178 double lon, lat;
179 bg::formula::sjoberg_intersection<double, bg::formula::andoyer_inverse, 1>
180 ::apply(-0.0872665, -0.0872665, -0.0872665, 0.0872665,
181 0.0, 1.57e-07, -0.392699, 1.57e-07,
182 lon, lat, bg::srs::spheroid<double>());
183 check_one("issue 612", lon, -0.087266500535674751);
184 check_one("issue 612", lat, 1.5892499139622920e-07);
185 }
186 }
187
test_main(int,char * [])188 int test_main(int, char*[])
189 {
190 for (size_t i = 0; i < expected_size; ++i)
191 {
192 test_all(expected[i]);
193 }
194
195 test_bugs();
196
197 return 0;
198 }
199