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