diff options
Diffstat (limited to 'boost/geometry/strategies/spherical/ssf.hpp')
-rw-r--r-- | boost/geometry/strategies/spherical/ssf.hpp | 109 |
1 files changed, 56 insertions, 53 deletions
diff --git a/boost/geometry/strategies/spherical/ssf.hpp b/boost/geometry/strategies/spherical/ssf.hpp index 41562950fd..81f3205e90 100644 --- a/boost/geometry/strategies/spherical/ssf.hpp +++ b/boost/geometry/strategies/spherical/ssf.hpp @@ -16,8 +16,9 @@ #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/radian_access.hpp> -#include <boost/geometry/util/select_coordinate_type.hpp> #include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/promote_floating_point.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> #include <boost/geometry/strategies/side.hpp> //#include <boost/geometry/strategies/concepts/side_concept.hpp> @@ -30,6 +31,43 @@ namespace boost { namespace geometry namespace strategy { namespace side { +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template <typename T> +int spherical_side_formula(T const& lambda1, T const& delta1, + T const& lambda2, T const& delta2, + T const& lambda, T const& delta) +{ + // Create temporary points (vectors) on unit a sphere + T const cos_delta1 = cos(delta1); + T const c1x = cos_delta1 * cos(lambda1); + T const c1y = cos_delta1 * sin(lambda1); + T const c1z = sin(delta1); + + T const cos_delta2 = cos(delta2); + T const c2x = cos_delta2 * cos(lambda2); + T const c2y = cos_delta2 * sin(lambda2); + T const c2z = sin(delta2); + + // (Third point is converted directly) + T const cos_delta = cos(delta); + + // Apply the "Spherical Side Formula" as presented on my blog + T const dist + = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda) + + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda) + + (c1x * c2y - c1y * c2x) * sin(delta); + + T zero = T(); + return dist > zero ? 1 + : dist < zero ? -1 + : 0; +} + +} +#endif // DOXYGEN_NO_DETAIL /*! \brief Check at which side of a Great Circle segment a point lies @@ -45,60 +83,25 @@ public : template <typename P1, typename P2, typename P> static inline int apply(P1 const& p1, P2 const& p2, P const& p) { - typedef typename boost::mpl::if_c + typedef typename promote_floating_point < - boost::is_void<CalculationType>::type::value, - - // Select at least a double... - typename select_most_precise + typename select_calculation_type_alt < - typename select_most_precise - < - typename select_most_precise - < - typename coordinate_type<P1>::type, - typename coordinate_type<P2>::type - >::type, - typename coordinate_type<P>::type - >::type, - double - >::type, - CalculationType - >::type coordinate_type; - - // Convenient shortcuts - typedef coordinate_type ct; - ct const lambda1 = get_as_radian<0>(p1); - ct const delta1 = get_as_radian<1>(p1); - ct const lambda2 = get_as_radian<0>(p2); - ct const delta2 = get_as_radian<1>(p2); - ct const lambda = get_as_radian<0>(p); - ct const delta = get_as_radian<1>(p); - - // Create temporary points (vectors) on unit a sphere - ct const cos_delta1 = cos(delta1); - ct const c1x = cos_delta1 * cos(lambda1); - ct const c1y = cos_delta1 * sin(lambda1); - ct const c1z = sin(delta1); - - ct const cos_delta2 = cos(delta2); - ct const c2x = cos_delta2 * cos(lambda2); - ct const c2y = cos_delta2 * sin(lambda2); - ct const c2z = sin(delta2); - - // (Third point is converted directly) - ct const cos_delta = cos(delta); - - // Apply the "Spherical Side Formula" as presented on my blog - ct const dist - = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda) - + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda) - + (c1x * c2y - c1y * c2x) * sin(delta); - - ct zero = ct(); - return dist > zero ? 1 - : dist < zero ? -1 - : 0; + CalculationType, + P1, P2, P + >::type + >::type calculation_type; + + calculation_type const lambda1 = get_as_radian<0>(p1); + calculation_type const delta1 = get_as_radian<1>(p1); + calculation_type const lambda2 = get_as_radian<0>(p2); + calculation_type const delta2 = get_as_radian<1>(p2); + calculation_type const lambda = get_as_radian<0>(p); + calculation_type const delta = get_as_radian<1>(p); + + return detail::spherical_side_formula(lambda1, delta1, + lambda2, delta2, + lambda, delta); } }; |