1 /*
2 * Copyright (C) 2010-2011 Dmitry Marakasov
3 *
4 * This file is part of glosm.
5 *
6 * glosm is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * glosm is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with glosm. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 #include <glosm/SphericalProjection.hh>
21
22 #include <glosm/geomath.h>
23
24 #include <cmath>
25
SphericalProjection()26 SphericalProjection::SphericalProjection() : Projection(&ProjectImpl, &UnProjectImpl) {
27 }
28
ProjectImpl(const Vector3i & point,const Vector3i & ref)29 Vector3f SphericalProjection::ProjectImpl(const Vector3i& point, const Vector3i& ref) {
30 double point_angle_x = ((double)point.x - (double)ref.x) * GEOM_DEG_TO_RAD;
31 double point_angle_y = (double)point.y * GEOM_DEG_TO_RAD;
32 double point_height = (double)point.z / GEOM_UNITSINMETER;
33
34 double ref_angle_y = (double)ref.y * GEOM_DEG_TO_RAD;
35 double ref_height = (double)ref.z / GEOM_UNITSINMETER;
36
37 #ifdef HAVE_SINCOS
38 double sinx, cosx, siny, cosy;
39 sincos(point_angle_x, &sinx, &cosx);
40 sincos(point_angle_y, &siny, &cosy);
41 #else
42 double cosx = cos(point_angle_x);
43 double cosy = cos(point_angle_y);
44 double sinx = sin(point_angle_x);
45 double siny = sin(point_angle_y);
46 #endif
47 Vector3d point_vector(
48 (WGS84_EARTH_EQ_RADIUS + point_height) * sinx * cosy,
49 (WGS84_EARTH_EQ_RADIUS + point_height) * siny,
50 (WGS84_EARTH_EQ_RADIUS + point_height) * cosx * cosy
51 );
52
53 #ifdef HAVE_SINCOS
54 double sinay, cosay;
55 sincos(ref_angle_y, &sinay, &cosay);
56 #else
57 double sinay = sin(ref_angle_y);
58 double cosay = cos(ref_angle_y);
59 #endif
60 return Vector3f(
61 point_vector.x,
62 point_vector.y * cosay - point_vector.z * sinay,
63 point_vector.y * sinay + point_vector.z * cosay - WGS84_EARTH_EQ_RADIUS - ref_height
64 );
65 }
66
UnProjectImpl(const Vector3f & point,const Vector3i & ref)67 Vector3i SphericalProjection::UnProjectImpl(const Vector3f& point, const Vector3i& ref) {
68 double ref_angle_y = (double)ref.y * GEOM_DEG_TO_RAD;
69 double ref_height = (double)ref.z / GEOM_UNITSINMETER;
70
71 #ifdef HAVE_SINCOS
72 double sinay, cosay;
73 sincos(ref_angle_y, &sinay, &cosay);
74 #else
75 double sinay = sin(ref_angle_y);
76 double cosay = cos(ref_angle_y);
77 #endif
78 Vector3d point_vector(
79 point.x,
80 sinay * (WGS84_EARTH_EQ_RADIUS + ref_height + point.z) + cosay * point.y,
81 cosay * (WGS84_EARTH_EQ_RADIUS + ref_height + point.z) - sinay * point.y
82 );
83
84 Vector3d spherical(
85 atan2(point_vector.x, point_vector.z),
86 atan2(point_vector.y, sqrt(point_vector.x * point_vector.x + point_vector.z * point_vector.z)),
87 point_vector.Length() - WGS84_EARTH_EQ_RADIUS
88 );
89
90 return Vector3i(
91 ref.x + (osmint_t)round(spherical.x * GEOM_RAD_TO_DEG),
92 (osmint_t)round(spherical.y * GEOM_RAD_TO_DEG),
93 (osmint_t)round(spherical.z * GEOM_UNITSINMETER)
94 );
95 }
96