summaryrefslogtreecommitdiff
path: root/boost/geometry/strategies/cartesian/side_by_triangle.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/strategies/cartesian/side_by_triangle.hpp')
-rw-r--r--boost/geometry/strategies/cartesian/side_by_triangle.hpp179
1 files changed, 157 insertions, 22 deletions
diff --git a/boost/geometry/strategies/cartesian/side_by_triangle.hpp b/boost/geometry/strategies/cartesian/side_by_triangle.hpp
index 5d589ffc86..77443d46a9 100644
--- a/boost/geometry/strategies/cartesian/side_by_triangle.hpp
+++ b/boost/geometry/strategies/cartesian/side_by_triangle.hpp
@@ -1,8 +1,13 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
-// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
-// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
+// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
+
+// This file was modified by Oracle on 2015.
+// Modifications copyright (c) 2015, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -22,6 +27,9 @@
#include <boost/geometry/util/select_coordinate_type.hpp>
#include <boost/geometry/strategies/side.hpp>
+#include <boost/geometry/algorithms/detail/relate/less.hpp>
+#include <boost/geometry/algorithms/detail/equals/point_point.hpp>
+
namespace boost { namespace geometry
{
@@ -38,6 +46,24 @@ namespace strategy { namespace side
template <typename CalculationType = void>
class side_by_triangle
{
+ template <typename Policy>
+ struct eps_policy
+ {
+ eps_policy() {}
+ template <typename Type>
+ eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
+ : policy(a, b, c, d)
+ {}
+ Policy policy;
+ };
+
+ struct eps_empty
+ {
+ eps_empty() {}
+ template <typename Type>
+ eps_empty(Type const&, Type const&, Type const&, Type const&) {}
+ };
+
public :
// Template member function, because it is not always trivial
@@ -47,23 +73,34 @@ public :
// Types can be all three different. Therefore it is
// not implemented (anymore) as "segment"
- template <typename coordinate_type, typename promoted_type, typename P1, typename P2, typename P>
- static inline promoted_type side_value(P1 const& p1, P2 const& p2, P const& p)
+ template
+ <
+ typename CoordinateType,
+ typename PromotedType,
+ typename P1,
+ typename P2,
+ typename P,
+ typename EpsPolicy
+ >
+ static inline
+ PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
{
- coordinate_type const x = get<0>(p);
- coordinate_type const y = get<1>(p);
+ CoordinateType const x = get<0>(p);
+ CoordinateType const y = get<1>(p);
- coordinate_type const sx1 = get<0>(p1);
- coordinate_type const sy1 = get<1>(p1);
- coordinate_type const sx2 = get<0>(p2);
- coordinate_type const sy2 = get<1>(p2);
+ CoordinateType const sx1 = get<0>(p1);
+ CoordinateType const sy1 = get<1>(p1);
+ CoordinateType const sx2 = get<0>(p2);
+ CoordinateType const sy2 = get<1>(p2);
- promoted_type const dx = sx2 - sx1;
- promoted_type const dy = sy2 - sy1;
- promoted_type const dpx = x - sx1;
- promoted_type const dpy = y - sy1;
+ PromotedType const dx = sx2 - sx1;
+ PromotedType const dy = sy2 - sy1;
+ PromotedType const dpx = x - sx1;
+ PromotedType const dpy = y - sy1;
- return geometry::detail::determinant<promoted_type>
+ eps_policy = EpsPolicy(dx, dy, dpx, dpy);
+
+ return geometry::detail::determinant<PromotedType>
(
dx, dy,
dpx, dpy
@@ -71,9 +108,99 @@ public :
}
+ template
+ <
+ typename CoordinateType,
+ typename PromotedType,
+ typename P1,
+ typename P2,
+ typename P
+ >
+ static inline
+ PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
+ {
+ eps_empty dummy;
+ return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
+ }
+
+
+ template
+ <
+ typename CoordinateType,
+ typename PromotedType,
+ bool AreAllIntegralCoordinates
+ >
+ struct compute_side_value
+ {
+ template <typename P1, typename P2, typename P, typename EpsPolicy>
+ static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
+ {
+ return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
+ }
+ };
+
+ template <typename CoordinateType, typename PromotedType>
+ struct compute_side_value<CoordinateType, PromotedType, false>
+ {
+ template <typename P1, typename P2, typename P, typename EpsPolicy>
+ static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
+ {
+ // For robustness purposes, first check if any two points are
+ // the same; in this case simply return that the points are
+ // collinear
+ if (geometry::detail::equals::equals_point_point(p1, p2)
+ || geometry::detail::equals::equals_point_point(p1, p)
+ || geometry::detail::equals::equals_point_point(p2, p))
+ {
+ return PromotedType(0);
+ }
+
+ // The side_by_triangle strategy computes the signed area of
+ // the point triplet (p1, p2, p); as such it is (in theory)
+ // invariant under cyclic permutations of its three arguments.
+ //
+ // In the context of numerical errors that arise in
+ // floating-point computations, and in order to make the strategy
+ // consistent with respect to cyclic permutations of its three
+ // arguments, we cyclically permute them so that the first
+ // argument is always the lexicographically smallest point.
+
+ geometry::detail::relate::less less;
+ if (less(p, p1))
+ {
+ if (less(p, p2))
+ {
+ // p is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
+ }
+ else
+ {
+ // p2 is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
+ }
+ }
+
+ if (less(p1, p2))
+ {
+ // p1 is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
+ }
+ else
+ {
+ // p2 is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
+ }
+ }
+ };
+
+
template <typename P1, typename P2, typename P>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
+ typedef typename coordinate_type<P1>::type coordinate_type1;
+ typedef typename coordinate_type<P2>::type coordinate_type2;
+ typedef typename coordinate_type<P>::type coordinate_type3;
+
typedef typename boost::mpl::if_c
<
boost::is_void<CalculationType>::type::value,
@@ -81,10 +208,9 @@ public :
<
typename select_most_precise
<
- typename coordinate_type<P1>::type,
- typename coordinate_type<P2>::type
+ coordinate_type1, coordinate_type2
>::type,
- typename coordinate_type<P>::type
+ coordinate_type3
>::type,
CalculationType
>::type coordinate_type;
@@ -96,10 +222,19 @@ public :
double
>::type promoted_type;
- promoted_type const s = side_value<coordinate_type, promoted_type>(p1, p2, p);
- promoted_type const zero = promoted_type();
+ bool const are_all_integral_coordinates =
+ boost::is_integral<coordinate_type1>::value
+ && boost::is_integral<coordinate_type2>::value
+ && boost::is_integral<coordinate_type3>::value;
- return math::equals(s, zero) ? 0
+ eps_policy< math::detail::equals_factor_policy<promoted_type> > epsp;
+ promoted_type s = compute_side_value
+ <
+ coordinate_type, promoted_type, are_all_integral_coordinates
+ >::apply(p1, p2, p, epsp);
+
+ promoted_type const zero = promoted_type();
+ return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
: s > zero ? 1
: -1;
}