summaryrefslogtreecommitdiff
path: root/boost/geometry/strategies/spherical/ssf.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/strategies/spherical/ssf.hpp')
-rw-r--r--boost/geometry/strategies/spherical/ssf.hpp109
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);
}
};