1 /* Copyright (c) 2017, Oracle and/or its affiliates. All rights reserved.
2
3 This program is free software; you can redistribute it and/or modify
4 it under the terms of the GNU General Public License, version 2.0,
5 as published by the Free Software Foundation.
6
7 This program is also distributed with certain software (including
8 but not limited to OpenSSL) that is licensed under separate terms,
9 as designated in a particular file or component or in included license
10 documentation. The authors of MySQL hereby grant you an additional
11 permission to link the program and your derivative works with the
12 separately licensed software that they have included with MySQL.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License, version 2.0, for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with this program; if not, write to the Free Software
21 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */
22
23 #ifndef UNITTEST_GUNIT_GIS_PROJECT_ON_POLE_H_INCLUDED
24 #define UNITTEST_GUNIT_GIS_PROJECT_ON_POLE_H_INCLUDED
25
26 #include <cmath> // M_PI_2
27
28 #include "sql/gis/geometries_cs.h" // gis::{Cartesian_*,Geographic_*}
29
30 namespace {
31
32 /**
33 * Convert to polar coordinates (-y as polar axis), interpret as sphere north
34 * pole azimuthal equidistant coordinates, then unproject to latiude /
35 * longtitude. This transformation does not preserve any properties, and in
36 * particular, straight lines do not become geodesics unless going through
37 * cartesian point (0,0), mapped to the north pole. Intersections should happen
38 * at cartesian (0,0) to be guaranteed preseved.
39 *
40 * A better projection would be something like polar gnomonic but for an
41 * ellipsoid surface. Such a projection would ensure straight lines get mapped
42 * to geodesics. However, such a projection seems impossible.
43 */
gis_project_on_pole(const gis::Cartesian_point & cart)44 gis::Geographic_point gis_project_on_pole(const gis::Cartesian_point &cart) {
45 using std::atan2;
46 using std::hypot;
47
48 auto x = cart.x();
49 auto y = cart.y();
50 // Input point should be inside radius pi/2 of origin. Otherwise
51 // clockwisedness of rings would not be preserved.
52 assert(hypot(x, y) <= M_PI_2);
53 // Rotate 45deg and project, so that -y points along the prime meridian and
54 // origin is on north pole, convert to polar, and translate origin to
55 // equator.
56 return {atan2(-y, x), (M_PI_2 - hypot(-y, x))};
57 }
58
gis_project_on_pole(const gis::Cartesian_linestring & cart)59 inline gis::Geographic_linestring gis_project_on_pole(
60 const gis::Cartesian_linestring &cart) {
61 gis::Geographic_linestring geom;
62 for (size_t i = 0; i < cart.size(); i++) {
63 geom.push_back(gis_project_on_pole(cart[i]));
64 }
65 return geom;
66 }
67
68 } // namespace
69
70 #endif // include guard
71