geometry/formulas/interpolate_point_spherical.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2019-2023 Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 #ifndef BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP
0014
0015 #include <boost/geometry/arithmetic/normalize.hpp>
0016 #include <boost/geometry/formulas/spherical.hpp>
0017 #include <boost/geometry/geometries/point.hpp>
0018
0019 namespace boost { namespace geometry { namespace formula
0020 {
0021
0022 template <typename CalculationType>
0023 class interpolate_point_spherical
0024 {
0025     typedef model::point<CalculationType, 3, cs::cartesian> point3d_t;
0026
0027 public :
0028
0029     template <typename Point>
0030     void compute_angle(Point const& p0,
0031                        Point const& p1,
0032                        CalculationType& angle01)
0033     {
0034         m_xyz0 = formula::sph_to_cart3d<point3d_t>(p0);
0035         m_xyz1 = formula::sph_to_cart3d<point3d_t>(p1);
0036         CalculationType const dot01 = geometry::dot_product(m_xyz0, m_xyz1);
0037         angle01 = acos(dot01);
0038     }
0039
0040     template <typename Point>
0041     void compute_axis(Point const& p0,
0042                       CalculationType const& angle01)
0043     {
0044         CalculationType const c0 = 0, c1 = 1;
0045         CalculationType const pi = math::pi<CalculationType>();
0046
0047         if (! math::equals(angle01, pi))
0048         {
0049             m_axis = geometry::cross_product(m_xyz0, m_xyz1);
0050             geometry::detail::vec_normalize(m_axis);
0051         }
0052         else // antipodal
0053         {
0054             CalculationType const half_pi = math::half_pi<CalculationType>();
0055             CalculationType const lat = geometry::get_as_radian<1>(p0);
0056
0057             if (math::equals(lat, half_pi))
0058             {
0059                 // pointing east, segment lies on prime meridian, going south
0060                 m_axis = point3d_t(c0, c1, c0);
0061             }
0062             else if (math::equals(lat, -half_pi))
0063             {
0064                 // pointing west, segment lies on prime meridian, going north
0065                 m_axis = point3d_t(c0, -c1, c0);
0066             }
0067             else
0068             {
0069                 // lon rotated west by pi/2 at equator
0070                 CalculationType const lon = geometry::get_as_radian<0>(p0);
0071                 m_axis = point3d_t(sin(lon), -cos(lon), c0);
0072             }
0073         }
0074     }
0075
0076     template <typename Point>
0077     void compute_point(CalculationType const& a, Point& p)
0078     {
0079         CalculationType const c1 = 1;
0080
0081         // Axis-Angle rotation
0082         // see: https://en.wikipedia.org/wiki/Axis-angle_representation
0083         CalculationType const cos_a = cos(a);
0084         CalculationType const sin_a = sin(a);
0085         // cos_a * v
0086         point3d_t s1 = m_xyz0;
0087         geometry::multiply_value(s1, cos_a);
0088         // sin_a * (n x v)
0089         point3d_t s2 = geometry::cross_product(m_axis, m_xyz0);
0090         geometry::multiply_value(s2, sin_a);
0091         // (1 - cos_a)(n.v) * n
0092         point3d_t s3 = m_axis;
0093         geometry::multiply_value(s3, (c1 - cos_a) *
0094                                  geometry::dot_product(m_axis, m_xyz0));
0095         // v_rot = cos_a * v + sin_a * (n x v) + (1 - cos_a)(n.v) * e
0096         point3d_t v_rot = s1;
0097         geometry::add_point(v_rot, s2);
0098         geometry::add_point(v_rot, s3);
0099
0100         p = formula::cart3d_to_sph<Point>(v_rot);
0101     }
0102
0103 private :
0104     point3d_t m_xyz0;
0105     point3d_t m_xyz1;
0106     point3d_t m_axis;
0107 };
0108
0109 }}} // namespace boost::geometry::formula
0110
0111 #endif // BOOST_GEOMETRY_FORMULAS_INTERPOLATE_POINT_SPHERICAL_HPP