boost/geometry/strategies/strategy_transform.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
   4 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
   5 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
   6
   7 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
   8 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
   9
  10 // Use, modification and distribution is subject to the Boost Software License,
  11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  12 // http://www.boost.org/LICENSE_1_0.txt)
  13
  14 #ifndef BOOST_GEOMETRY_STRATEGIES_STRATEGY_TRANSFORM_HPP
  15 #define BOOST_GEOMETRY_STRATEGIES_STRATEGY_TRANSFORM_HPP
  16
  17 #include <cstddef>
  18 #include <cmath>
  19 #include <functional>
  20
  21 #include <boost/numeric/conversion/cast.hpp>
  22
  23 #include <boost/geometry/algorithms/convert.hpp>
  24 #include <boost/geometry/arithmetic/arithmetic.hpp>
  25 #include <boost/geometry/core/access.hpp>
  26 #include <boost/geometry/core/radian_access.hpp>
  27 #include <boost/geometry/core/coordinate_dimension.hpp>
  28 #include <boost/geometry/strategies/transform.hpp>
  29
  30 #include <boost/geometry/util/math.hpp>
  31 #include <boost/geometry/util/select_coordinate_type.hpp>
  32
  33 namespace boost { namespace geometry
  34 {
  35
  36 namespace strategy { namespace transform
  37 {
  38
  39 #ifndef DOXYGEN_NO_DETAIL
  40 namespace detail
  41 {
  42
  43 template
  44 <
  45     typename Src, typename Dst,
  46     std::size_t D, std::size_t N,
  47     template <typename> class F
  48 >
  49 struct transform_coordinates
  50 {
  51     template <typename T>
  52     static inline void transform(Src const& source, Dst& dest, T value)
  53     {
  54         typedef typename select_coordinate_type<Src, Dst>::type coordinate_type;
  55
  56         F<coordinate_type> function;
  57         set<D>(dest, boost::numeric_cast<coordinate_type>(function(get<D>(source), value)));
  58         transform_coordinates<Src, Dst, D + 1, N, F>::transform(source, dest, value);
  59     }
  60 };
  61
  62 template
  63 <
  64     typename Src, typename Dst,
  65     std::size_t N,
  66     template <typename> class F
  67 >
  68 struct transform_coordinates<Src, Dst, N, N, F>
  69 {
  70     template <typename T>
  71     static inline void transform(Src const& , Dst& , T )
  72     {
  73     }
  74 };
  75
  76 } // namespace detail
  77 #endif // DOXYGEN_NO_DETAIL
  78
  79
  80 /*!
  81     \brief Transformation strategy to copy one point to another using assignment operator
  82     \ingroup transform
  83     \tparam P point type
  84  */
  85 template <typename P>
  86 struct copy_direct
  87 {
  88     inline bool apply(P const& p1, P& p2) const
  89     {
  90         p2 = p1;
  91         return true;
  92     }
  93 };
  94
  95 /*!
  96     \brief Transformation strategy to do copy a point, copying per coordinate.
  97     \ingroup transform
  98     \tparam P1 first point type
  99     \tparam P2 second point type
 100  */
 101 template <typename P1, typename P2>
 102 struct copy_per_coordinate
 103 {
 104     inline bool apply(P1 const& p1, P2& p2) const
 105     {
 106         // Defensive check, dimensions are equal, selected by specialization
 107         assert_dimension_equal<P1, P2>();
 108
 109         geometry::convert(p1, p2);
 110         return true;
 111     }
 112 };
 113
 114
 115 /*!
 116     \brief Transformation strategy to go from degree to radian and back
 117     \ingroup transform
 118     \tparam P1 first point type
 119     \tparam P2 second point type
 120     \tparam F additional functor to divide or multiply with d2r
 121  */
 122 template <typename P1, typename P2, template <typename> class F>
 123 struct degree_radian_vv
 124 {
 125     inline bool apply(P1 const& p1, P2& p2) const
 126     {
 127         // Spherical coordinates always have 2 coordinates measured in angles
 128         // The optional third one is distance/height, provided in another strategy
 129         // Polar coordinates having one angle, will be also in another strategy
 130         assert_dimension<P1, 2>();
 131         assert_dimension<P2, 2>();
 132
 133         detail::transform_coordinates<P1, P2, 0, 2, F>::transform(p1, p2, math::d2r);
 134         return true;
 135     }
 136 };
 137
 138 template <typename P1, typename P2, template <typename> class F>
 139 struct degree_radian_vv_3
 140 {
 141     inline bool apply(P1 const& p1, P2& p2) const
 142     {
 143         assert_dimension<P1, 3>();
 144         assert_dimension<P2, 3>();
 145
 146         detail::transform_coordinates<P1, P2, 0, 2, F>::transform(p1, p2, math::d2r);
 147         // Copy height or other third dimension
 148         set<2>(p2, get<2>(p1));
 149         return true;
 150     }
 151 };
 152
 153
 154 #ifndef DOXYGEN_NO_DETAIL
 155 namespace detail
 156 {
 157
 158     /// Helper function for conversion, phi/theta are in radians
 159     template <typename P, typename T, typename R>
 160     inline void spherical_polar_to_cartesian(T phi, T theta, R r, P& p)
 161     {
 162         assert_dimension<P, 3>();
 163
 164         // http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_spherical_coordinates
 165         // http://www.vias.org/comp_geometry/math_coord_convert_3d.htm
 166         // https://moodle.polymtl.ca/file.php/1183/Autres_Documents/Derivation_for_Spherical_Co-ordinates.pdf
 167         // http://en.citizendium.org/wiki/Spherical_polar_coordinates
 168
 169         // Phi = first, theta is second, r is third, see documentation on cs::spherical
 170
 171         // (calculations are splitted to implement ttmath)
 172
 173         T r_sin_theta = r;
 174         T r_cos_theta = r;
 175         r_sin_theta *= sin(theta);
 176         r_cos_theta *= cos(theta);
 177
 178         set<0>(p, r_sin_theta * cos(phi));
 179         set<1>(p, r_sin_theta * sin(phi));
 180         set<2>(p, r_cos_theta);
 181     }
 182
 183     /// Helper function for conversion, lambda/delta (lon lat) are in radians
 184     template <typename P, typename T, typename R>
 185     inline void spherical_equatorial_to_cartesian(T lambda, T delta, R r, P& p)
 186     {
 187         assert_dimension<P, 3>();
 188
 189         // http://mathworld.wolfram.com/GreatCircle.html
 190         // http://www.spenvis.oma.be/help/background/coortran/coortran.html WRONG
 191
 192         T r_cos_delta = r;
 193         T r_sin_delta = r;
 194         r_cos_delta *= cos(delta);
 195         r_sin_delta *= sin(delta);
 196
 197         set<0>(p, r_cos_delta * cos(lambda));
 198         set<1>(p, r_cos_delta * sin(lambda));
 199         set<2>(p, r_sin_delta);
 200     }
 201
 202
 203     /// Helper function for conversion
 204     template <typename P, typename T>
 205     inline bool cartesian_to_spherical2(T x, T y, T z, P& p)
 206     {
 207         assert_dimension<P, 2>();
 208
 209         // http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_Cartesian_coordinates
 210
 211 #if defined(BOOST_GEOMETRY_TRANSFORM_CHECK_UNIT_SPHERE)
 212         // TODO: MAYBE ONLY IF TO BE CHECKED?
 213         T const r = /*sqrt not necessary, sqrt(1)=1*/ (x * x + y * y + z * z);
 214
 215         // Unit sphere, so r should be 1
 216         if (geometry::math::abs(r - 1.0) > T(1e-6))
 217         {
 218             return false;
 219         }
 220         // end todo
 221 #endif
 222
 223         set_from_radian<0>(p, atan2(y, x));
 224         set_from_radian<1>(p, acos(z));
 225         return true;
 226     }
 227
 228     template <typename P, typename T>
 229     inline bool cartesian_to_spherical_equatorial2(T x, T y, T z, P& p)
 230     {
 231         assert_dimension<P, 2>();
 232
 233         set_from_radian<0>(p, atan2(y, x));
 234         set_from_radian<1>(p, asin(z));
 235         return true;
 236     }
 237
 238
 239     template <typename P, typename T>
 240     inline bool cartesian_to_spherical3(T x, T y, T z, P& p)
 241     {
 242         assert_dimension<P, 3>();
 243
 244         // http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_Cartesian_coordinates
 245         T const r = math::sqrt(x * x + y * y + z * z);
 246         set<2>(p, r);
 247         set_from_radian<0>(p, atan2(y, x));
 248         if (r > 0.0)
 249         {
 250             set_from_radian<1>(p, acos(z / r));
 251             return true;
 252         }
 253         return false;
 254     }
 255
 256     template <typename P, typename T>
 257     inline bool cartesian_to_spherical_equatorial3(T x, T y, T z, P& p)
 258     {
 259         assert_dimension<P, 3>();
 260
 261         // http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_Cartesian_coordinates
 262         T const r = math::sqrt(x * x + y * y + z * z);
 263         set<2>(p, r);
 264         set_from_radian<0>(p, atan2(y, x));
 265         if (r > 0.0)
 266         {
 267             set_from_radian<1>(p, asin(z / r));
 268             return true;
 269         }
 270         return false;
 271     }
 272
 273 } // namespace detail
 274 #endif // DOXYGEN_NO_DETAIL
 275
 276
 277 /*!
 278     \brief Transformation strategy for 2D spherical (phi,theta) to 3D cartesian (x,y,z)
 279     \details on Unit sphere
 280     \ingroup transform
 281     \tparam P1 first point type
 282     \tparam P2 second point type
 283  */
 284 template <typename P1, typename P2>
 285 struct from_spherical_polar_2_to_cartesian_3
 286 {
 287     inline bool apply(P1 const& p1, P2& p2) const
 288     {
 289         assert_dimension<P1, 2>();
 290         detail::spherical_polar_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2);
 291         return true;
 292     }
 293 };
 294
 295 template <typename P1, typename P2>
 296 struct from_spherical_equatorial_2_to_cartesian_3
 297 {
 298     inline bool apply(P1 const& p1, P2& p2) const
 299     {
 300         assert_dimension<P1, 2>();
 301         detail::spherical_equatorial_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2);
 302         return true;
 303     }
 304 };
 305
 306
 307 /*!
 308     \brief Transformation strategy for 3D spherical (phi,theta,r) to 3D cartesian (x,y,z)
 309     \ingroup transform
 310     \tparam P1 first point type
 311     \tparam P2 second point type
 312  */
 313 template <typename P1, typename P2>
 314 struct from_spherical_polar_3_to_cartesian_3
 315 {
 316     inline bool apply(P1 const& p1, P2& p2) const
 317     {
 318         assert_dimension<P1, 3>();
 319         detail::spherical_polar_to_cartesian(
 320                     get_as_radian<0>(p1), get_as_radian<1>(p1), get<2>(p1), p2);
 321         return true;
 322     }
 323 };
 324
 325 template <typename P1, typename P2>
 326 struct from_spherical_equatorial_3_to_cartesian_3
 327 {
 328     inline bool apply(P1 const& p1, P2& p2) const
 329     {
 330         assert_dimension<P1, 3>();
 331         detail::spherical_equatorial_to_cartesian(
 332                     get_as_radian<0>(p1), get_as_radian<1>(p1), get<2>(p1), p2);
 333         return true;
 334     }
 335 };
 336
 337
 338 /*!
 339     \brief Transformation strategy for 3D cartesian (x,y,z) to 2D spherical (phi,theta)
 340     \details on Unit sphere
 341     \ingroup transform
 342     \tparam P1 first point type
 343     \tparam P2 second point type
 344     \note If x,y,z point is not lying on unit sphere, transformation will return false
 345  */
 346 template <typename P1, typename P2>
 347 struct from_cartesian_3_to_spherical_polar_2
 348 {
 349     inline bool apply(P1 const& p1, P2& p2) const
 350     {
 351         assert_dimension<P1, 3>();
 352         return detail::cartesian_to_spherical2(get<0>(p1), get<1>(p1), get<2>(p1), p2);
 353     }
 354 };
 355
 356 template <typename P1, typename P2>
 357 struct from_cartesian_3_to_spherical_equatorial_2
 358 {
 359     inline bool apply(P1 const& p1, P2& p2) const
 360     {
 361         assert_dimension<P1, 3>();
 362         return detail::cartesian_to_spherical_equatorial2(get<0>(p1), get<1>(p1), get<2>(p1), p2);
 363     }
 364 };
 365
 366
 367 /*!
 368     \brief Transformation strategy for 3D cartesian (x,y,z) to 3D spherical (phi,theta,r)
 369     \ingroup transform
 370     \tparam P1 first point type
 371     \tparam P2 second point type
 372  */
 373 template <typename P1, typename P2>
 374 struct from_cartesian_3_to_spherical_polar_3
 375 {
 376     inline bool apply(P1 const& p1, P2& p2) const
 377     {
 378         assert_dimension<P1, 3>();
 379         return detail::cartesian_to_spherical3(get<0>(p1), get<1>(p1), get<2>(p1), p2);
 380     }
 381 };
 382
 383 template <typename P1, typename P2>
 384 struct from_cartesian_3_to_spherical_equatorial_3
 385 {
 386     inline bool apply(P1 const& p1, P2& p2) const
 387     {
 388         assert_dimension<P1, 3>();
 389         return detail::cartesian_to_spherical_equatorial3(get<0>(p1), get<1>(p1), get<2>(p1), p2);
 390     }
 391 };
 392
 393 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 394
 395 namespace services
 396 {
 397
 398 /// Specialization for same coordinate system family, same system, same dimension, same point type, can be copied
 399 template <typename CoordSysTag, typename CoordSys, std::size_t D, typename P>
 400 struct default_strategy<CoordSysTag, CoordSysTag, CoordSys, CoordSys, D, D, P, P>
 401 {
 402     typedef copy_direct<P> type;
 403 };
 404
 405 /// Specialization for same coordinate system family and system, same dimension, different point type, copy per coordinate
 406 template <typename CoordSysTag, typename CoordSys, std::size_t D, typename P1, typename P2>
 407 struct default_strategy<CoordSysTag, CoordSysTag, CoordSys, CoordSys, D, D, P1, P2>
 408 {
 409     typedef copy_per_coordinate<P1, P2> type;
 410 };
 411
 412 /// Specialization to transform from degree to radian for any coordinate system / point type combination
 413 template <typename CoordSysTag, template<typename> class CoordSys, typename P1, typename P2>
 414 struct default_strategy<CoordSysTag, CoordSysTag, CoordSys<degree>, CoordSys<radian>, 2, 2, P1, P2>
 415 {
 416     typedef degree_radian_vv<P1, P2, std::multiplies> type;
 417 };
 418
 419 /// Specialization to transform from radian to degree for any coordinate system / point type combination
 420 template <typename CoordSysTag, template<typename> class CoordSys, typename P1, typename P2>
 421 struct default_strategy<CoordSysTag, CoordSysTag, CoordSys<radian>, CoordSys<degree>, 2, 2, P1, P2>
 422 {
 423     typedef degree_radian_vv<P1, P2, std::divides> type;
 424 };
 425
 426
 427 /// Specialization degree->radian in 3D
 428 template <typename CoordSysTag, template<typename> class CoordSys, typename P1, typename P2>
 429 struct default_strategy<CoordSysTag, CoordSysTag, CoordSys<degree>, CoordSys<radian>, 3, 3, P1, P2>
 430 {
 431     typedef degree_radian_vv_3<P1, P2, std::multiplies> type;
 432 };
 433
 434 /// Specialization radian->degree in 3D
 435 template <typename CoordSysTag, template<typename> class CoordSys, typename P1, typename P2>
 436 struct default_strategy<CoordSysTag, CoordSysTag, CoordSys<radian>, CoordSys<degree>, 3, 3, P1, P2>
 437 {
 438     typedef degree_radian_vv_3<P1, P2, std::divides> type;
 439 };
 440
 441 /// Specialization to transform from unit sphere(phi,theta) to XYZ
 442 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 443 struct default_strategy<spherical_polar_tag, cartesian_tag, CoordSys1, CoordSys2, 2, 3, P1, P2>
 444 {
 445     typedef from_spherical_polar_2_to_cartesian_3<P1, P2> type;
 446 };
 447
 448 /// Specialization to transform from sphere(phi,theta,r) to XYZ
 449 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 450 struct default_strategy<spherical_polar_tag, cartesian_tag, CoordSys1, CoordSys2, 3, 3, P1, P2>
 451 {
 452     typedef from_spherical_polar_3_to_cartesian_3<P1, P2> type;
 453 };
 454
 455 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 456 struct default_strategy<spherical_equatorial_tag, cartesian_tag, CoordSys1, CoordSys2, 2, 3, P1, P2>
 457 {
 458     typedef from_spherical_equatorial_2_to_cartesian_3<P1, P2> type;
 459 };
 460
 461 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 462 struct default_strategy<spherical_equatorial_tag, cartesian_tag, CoordSys1, CoordSys2, 3, 3, P1, P2>
 463 {
 464     typedef from_spherical_equatorial_3_to_cartesian_3<P1, P2> type;
 465 };
 466
 467 /// Specialization to transform from XYZ to unit sphere(phi,theta)
 468 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 469 struct default_strategy<cartesian_tag, spherical_polar_tag, CoordSys1, CoordSys2, 3, 2, P1, P2>
 470 {
 471     typedef from_cartesian_3_to_spherical_polar_2<P1, P2> type;
 472 };
 473
 474 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 475 struct default_strategy<cartesian_tag, spherical_equatorial_tag, CoordSys1, CoordSys2, 3, 2, P1, P2>
 476 {
 477     typedef from_cartesian_3_to_spherical_equatorial_2<P1, P2> type;
 478 };
 479
 480 /// Specialization to transform from XYZ to sphere(phi,theta,r)
 481 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 482 struct default_strategy<cartesian_tag, spherical_polar_tag, CoordSys1, CoordSys2, 3, 3, P1, P2>
 483 {
 484     typedef from_cartesian_3_to_spherical_polar_3<P1, P2> type;
 485 };
 486 template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
 487 struct default_strategy<cartesian_tag, spherical_equatorial_tag, CoordSys1, CoordSys2, 3, 3, P1, P2>
 488 {
 489     typedef from_cartesian_3_to_spherical_equatorial_3<P1, P2> type;
 490 };
 491
 492
 493 } // namespace services
 494
 495
 496 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 497
 498
 499 }} // namespace strategy::transform
 500
 501
 502 }} // namespace boost::geometry
 503
 504 #endif // BOOST_GEOMETRY_STRATEGIES_STRATEGY_TRANSFORM_HPP