1 // Boost.Geometry
2 
3 // Copyright (c) 2019 Oracle and/or its affiliates.
4 
5 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
6 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
7 
8 // Use, modification and distribution is subject to the Boost Software License,
9 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
10 // http://www.boost.org/LICENSE_1_0.txt)
11 
12 #ifndef BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
13 #define BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
14 
15 namespace boost { namespace geometry { namespace formula
16 {
17 
18 template <typename CalculationType>
19 class interpolate_point_spherical
20 {
21     typedef model::point<CalculationType, 3, cs::cartesian> point3d_t;
22 
23 public :
24 
25     template <typename Point>
compute_angle(Point const & p0,Point const & p1,CalculationType & angle01)26     void compute_angle(Point const& p0,
27                        Point const& p1,
28                        CalculationType& angle01)
29     {
30         m_xyz0 = formula::sph_to_cart3d<point3d_t>(p0);
31         m_xyz1 = formula::sph_to_cart3d<point3d_t>(p1);
32         CalculationType const dot01 = geometry::dot_product(m_xyz0, m_xyz1);
33         angle01 = acos(dot01);
34     }
35 
36     template <typename Point>
compute_axis(Point const & p0,CalculationType const & angle01)37     void compute_axis(Point const& p0,
38                       CalculationType const& angle01)
39     {
40         CalculationType const c0 = 0, c1 = 1;
41         CalculationType const pi = math::pi<CalculationType>();
42 
43         if (! math::equals(angle01, pi))
44         {
45             m_axis = geometry::cross_product(m_xyz0, m_xyz1);
46             geometry::detail::vec_normalize(m_axis);
47         }
48         else // antipodal
49         {
50             CalculationType const half_pi = math::half_pi<CalculationType>();
51             CalculationType const lat = geometry::get_as_radian<1>(p0);
52 
53             if (math::equals(lat, half_pi))
54             {
55                 // pointing east, segment lies on prime meridian, going south
56                 m_axis = point3d_t(c0, c1, c0);
57             }
58             else if (math::equals(lat, -half_pi))
59             {
60                 // pointing west, segment lies on prime meridian, going north
61                 m_axis = point3d_t(c0, -c1, c0);
62             }
63             else
64             {
65                 // lon rotated west by pi/2 at equator
66                 CalculationType const lon = geometry::get_as_radian<0>(p0);
67                 m_axis = point3d_t(sin(lon), -cos(lon), c0);
68             }
69         }
70     }
71 
72     template <typename Point>
compute_point(CalculationType const & a,Point & p)73     void compute_point(CalculationType const& a, Point& p)
74     {
75         CalculationType const c1 = 1;
76 
77         // Axis-Angle rotation
78         // see: https://en.wikipedia.org/wiki/Axis-angle_representation
79         CalculationType const cos_a = cos(a);
80         CalculationType const sin_a = sin(a);
81         // cos_a * v
82         point3d_t s1 = m_xyz0;
83         geometry::multiply_value(s1, cos_a);
84         // sin_a * (n x v)
85         point3d_t s2 = geometry::cross_product(m_axis, m_xyz0);
86         geometry::multiply_value(s2, sin_a);
87         // (1 - cos_a)(n.v) * n
88         point3d_t s3 = m_axis;
89         geometry::multiply_value(s3, (c1 - cos_a) *
90                                  geometry::dot_product(m_axis, m_xyz0));
91         // v_rot = cos_a * v + sin_a * (n x v) + (1 - cos_a)(n.v) * e
92         point3d_t v_rot = s1;
93         geometry::add_point(v_rot, s2);
94         geometry::add_point(v_rot, s3);
95 
96         p = formula::cart3d_to_sph<Point>(v_rot);
97     }
98 
99 private :
100     point3d_t m_xyz0;
101     point3d_t m_xyz1;
102     point3d_t m_axis;
103 };
104 
105 }}} // namespace boost::geometry::formula
106 
107 #endif // BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
108