diff options
Diffstat (limited to 'boost/geometry/strategies/cartesian')
9 files changed, 491 insertions, 79 deletions
diff --git a/boost/geometry/strategies/cartesian/buffer_end_round.hpp b/boost/geometry/strategies/cartesian/buffer_end_round.hpp index 74780d6165..a233f1c4be 100644 --- a/boost/geometry/strategies/cartesian/buffer_end_round.hpp +++ b/boost/geometry/strategies/cartesian/buffer_end_round.hpp @@ -95,8 +95,9 @@ public : //! \brief Constructs the strategy //! \param points_per_circle points which would be used for a full circle + //! (if points_per_circle is smaller than 4, it is internally set to 4) explicit inline end_round(std::size_t points_per_circle = 90) - : m_points_per_circle(points_per_circle) + : m_points_per_circle((points_per_circle < 4u) ? 4u : points_per_circle) {} #ifndef DOXYGEN_SHOULD_SKIP_THIS @@ -106,7 +107,7 @@ public : inline void apply(Point const& penultimate_point, Point const& perp_left_point, Point const& ultimate_point, - Point const& , + Point const& perp_right_point, buffer_side_selector side, DistanceStrategy const& distance, RangeOut& range_out) const @@ -142,6 +143,13 @@ public : set<1>(shifted_point, get<1>(ultimate_point) + dist_half_diff * sin(alpha)); generate_points(shifted_point, alpha, (dist_left + dist_right) / two, range_out); } + + if (m_points_per_circle % 2 == 1) + { + // For a half circle, if the number of points is not even, + // we should insert the end point too, to generate a full cap + range_out.push_back(perp_right_point); + } } template <typename NumericType> diff --git a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp index 8fcf3b996c..99ec80527f 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp @@ -35,6 +35,8 @@ namespace strategy { namespace buffer their length. The miter is not changed to a bevel form (as done in some other software), it is just adapted to the specified miter_limit but keeps its miter form. + If the buffer distance is 5.0, and the miter limit is 2.0, generated points + will be located at a distance of at most 10.0 (2*5) units. This strategy is only applicable for Cartesian coordinate systems. \qbk{ diff --git a/boost/geometry/strategies/cartesian/buffer_join_round.hpp b/boost/geometry/strategies/cartesian/buffer_join_round.hpp index 9e467c85a0..9ec51cd1ec 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_round.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_round.hpp @@ -1,6 +1,6 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2012-2015 Barend Gehrels, Amsterdam, the Netherlands. // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -9,6 +9,8 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_HPP +#include <algorithm> + #include <boost/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/policies/compare.hpp> @@ -69,34 +71,38 @@ private : DistanceType const& buffer_distance, RangeOut& range_out) const { - PromotedType dx1 = get<0>(perp1) - get<0>(vertex); - PromotedType dy1 = get<1>(perp1) - get<1>(vertex); - PromotedType dx2 = get<0>(perp2) - get<0>(vertex); - PromotedType dy2 = get<1>(perp2) - get<1>(vertex); + PromotedType const dx1 = get<0>(perp1) - get<0>(vertex); + PromotedType const dy1 = get<1>(perp1) - get<1>(vertex); + PromotedType const dx2 = get<0>(perp2) - get<0>(vertex); + PromotedType const dy2 = get<1>(perp2) - get<1>(vertex); - BOOST_ASSERT(buffer_distance != 0); + PromotedType const two = 2.0; + PromotedType const two_pi = two * geometry::math::pi<PromotedType>(); - dx1 /= buffer_distance; - dy1 /= buffer_distance; - dx2 /= buffer_distance; - dy2 /= buffer_distance; + PromotedType const angle1 = atan2(dy1, dx1); + PromotedType angle2 = atan2(dy2, dx2); + while (angle2 > angle1) + { + angle2 -= two_pi; + } + PromotedType const angle_diff = angle1 - angle2; - PromotedType angle_diff = acos(dx1 * dx2 + dy1 * dy2); + // Divide the angle into an integer amount of steps to make it + // visually correct also for a low number of points / circle - PromotedType two = 2.0; - PromotedType steps = m_points_per_circle; - int n = boost::numeric_cast<int>(steps * angle_diff - / (two * geometry::math::pi<PromotedType>())); + // If a full circle is divided into 3 parts (e.g. angle is 125), + // the one point in between must still be generated + // The calculation below: + // - generates 1 point in between for an angle of 125 based on 3 points + // - generates 0 points in between for an angle of 90 based on 4 points - if (n <= 1) - { - return; - } + int const n = (std::max)(static_cast<int>( + ceil(m_points_per_circle * angle_diff / two_pi)), 1); - PromotedType const angle1 = atan2(dy1, dx1); - PromotedType diff = angle_diff / PromotedType(n); + PromotedType const diff = angle_diff / static_cast<PromotedType>(n); PromotedType a = angle1 - diff; + // Walk to n - 1 to avoid generating the last point for (int i = 0; i < n - 1; i++, a -= diff) { Point p; diff --git a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp index f64a82d8fc..86ebc43c9c 100644 --- a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp +++ b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp @@ -1,5 +1,12 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands. + +// Copyright (c) 2012-2015 Barend Gehrels, Amsterdam, the Netherlands. + +// 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 + // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) @@ -45,9 +52,10 @@ class point_circle { public : //! \brief Constructs the strategy - //! \param count number of points for the created circle + //! \param count number of points for the created circle (if count + //! is smaller than 3, count is internally set to 3) explicit point_circle(std::size_t count = 90) - : m_count(count) + : m_count((count < 3u) ? 3u : count) {} #ifndef DOXYGEN_SHOULD_SKIP_THIS diff --git a/boost/geometry/strategies/cartesian/cart_intersect.hpp b/boost/geometry/strategies/cartesian/cart_intersect.hpp index 66af2d2e9c..a7bd385226 100644 --- a/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -1,7 +1,13 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. +// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland. + +// This file was modified by Oracle on 2014. +// Modifications copyright (c) 2014, Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -194,7 +200,11 @@ struct relate_cartesian_segments get<1>(robust_b1) - get<1>(robust_a1), robust_db0, robust_db); - if (robust_da0 == 0) + math::detail::equals_factor_policy<robust_coordinate_type> + policy(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b); + robust_coordinate_type const zero = 0; + if (math::detail::equals_by_policy(robust_da0, zero, policy) + || math::detail::equals_by_policy(robust_db0, zero, policy)) { // If this is the case, no rescaling is done for FP precision. // We set it to collinear, but it indicates a robustness issue. @@ -211,25 +221,31 @@ struct relate_cartesian_segments if (collinear) { - bool const collinear_use_first - = geometry::math::abs(robust_dx_a) + geometry::math::abs(robust_dx_b) - >= geometry::math::abs(robust_dy_a) + geometry::math::abs(robust_dy_b); - - // Degenerate cases: segments of single point, lying on other segment, are not disjoint - // This situation is collinear too - - if (collinear_use_first) + std::pair<bool, bool> const collinear_use_first + = is_x_more_significant(geometry::math::abs(robust_dx_a), + geometry::math::abs(robust_dy_a), + geometry::math::abs(robust_dx_b), + geometry::math::abs(robust_dy_b), + a_is_point, b_is_point); + + if (collinear_use_first.second) { - return relate_collinear<0, ratio_type>(a, b, - robust_a1, robust_a2, robust_b1, robust_b2, - a_is_point, b_is_point); - } - else - { - // Y direction contains larger segments (maybe dx is zero) - return relate_collinear<1, ratio_type>(a, b, - robust_a1, robust_a2, robust_b1, robust_b2, - a_is_point, b_is_point); + // Degenerate cases: segments of single point, lying on other segment, are not disjoint + // This situation is collinear too + + if (collinear_use_first.first) + { + return relate_collinear<0, ratio_type>(a, b, + robust_a1, robust_a2, robust_b1, robust_b2, + a_is_point, b_is_point); + } + else + { + // Y direction contains larger segments (maybe dx is zero) + return relate_collinear<1, ratio_type>(a, b, + robust_a1, robust_a2, robust_b1, robust_b2, + a_is_point, b_is_point); + } } } @@ -237,6 +253,40 @@ struct relate_cartesian_segments } private: + // first is true if x is more significant + // second is true if the more significant difference is not 0 + template <typename RobustCoordinateType> + static inline std::pair<bool, bool> + is_x_more_significant(RobustCoordinateType const& abs_robust_dx_a, + RobustCoordinateType const& abs_robust_dy_a, + RobustCoordinateType const& abs_robust_dx_b, + RobustCoordinateType const& abs_robust_dy_b, + bool const a_is_point, + bool const b_is_point) + { + //BOOST_ASSERT_MSG(!(a_is_point && b_is_point), "both segments shouldn't be degenerated"); + + // for degenerated segments the second is always true because this function + // shouldn't be called if both segments were degenerated + + if (a_is_point) + { + return std::make_pair(abs_robust_dx_b >= abs_robust_dy_b, true); + } + else if (b_is_point) + { + return std::make_pair(abs_robust_dx_a >= abs_robust_dy_a, true); + } + else + { + RobustCoordinateType const min_dx = (std::min)(abs_robust_dx_a, abs_robust_dx_b); + RobustCoordinateType const min_dy = (std::min)(abs_robust_dy_a, abs_robust_dy_b); + return min_dx == min_dy ? + std::make_pair(true, min_dx > RobustCoordinateType(0)) : + std::make_pair(min_dx > min_dy, true); + } + } + template < std::size_t Dimension, @@ -319,17 +369,58 @@ private: RobustType const length_a = oa_2 - oa_1; // no abs, see above RobustType const length_b = ob_2 - ob_1; - RatioType const ra_from(oa_1 - ob_1, length_b); - RatioType const ra_to(oa_2 - ob_1, length_b); - RatioType const rb_from(ob_1 - oa_1, length_a); - RatioType const rb_to(ob_2 - oa_1, length_a); + RatioType ra_from(oa_1 - ob_1, length_b); + RatioType ra_to(oa_2 - ob_1, length_b); + RatioType rb_from(ob_1 - oa_1, length_a); + RatioType rb_to(ob_2 - oa_1, length_a); + + // use absolute measure to detect endpoints intersection + // NOTE: it'd be possible to calculate bx_wrt_a using ax_wrt_b values + int const a1_wrt_b = position_value(oa_1, ob_1, ob_2); + int const a2_wrt_b = position_value(oa_2, ob_1, ob_2); + int const b1_wrt_a = position_value(ob_1, oa_1, oa_2); + int const b2_wrt_a = position_value(ob_2, oa_1, oa_2); + + // fix the ratios if necessary + // CONSIDER: fixing ratios also in other cases, if they're inconsistent + // e.g. if ratio == 1 or 0 (so IP at the endpoint) + // but position value indicates that the IP is in the middle of the segment + // because one of the segments is very long + // In such case the ratios could be moved into the middle direction + // by some small value (e.g. EPS+1ULP) + if (a1_wrt_b == 1) + { + ra_from.assign(0, 1); + rb_from.assign(0, 1); + } + else if (a1_wrt_b == 3) + { + ra_from.assign(1, 1); + rb_to.assign(0, 1); + } + + if (a2_wrt_b == 1) + { + ra_to.assign(0, 1); + rb_from.assign(1, 1); + } + else if (a2_wrt_b == 3) + { + ra_to.assign(1, 1); + rb_to.assign(1, 1); + } - if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right())) + if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3)) + //if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right())) { return Policy::disjoint(); } - return Policy::segments_collinear(a, b, ra_from, ra_to, rb_from, rb_to); + bool const opposite = math::sign(length_a) != math::sign(length_b); + + return Policy::segments_collinear(a, b, opposite, + a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a, + ra_from, ra_to, rb_from, rb_to); } /// Relate segments where one is degenerate @@ -351,8 +442,32 @@ private: // b1/b2 (4..4) // Ratio: (4-2)/(6-2) RatioType const ratio(d - s1, s2 - s1); + + if (!ratio.on_segment()) + { + return Policy::disjoint(); + } + return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate); } + + template <typename ProjCoord1, typename ProjCoord2> + static inline int position_value(ProjCoord1 const& ca1, + ProjCoord2 const& cb1, + ProjCoord2 const& cb2) + { + // S1x 0 1 2 3 4 + // S2 |----------> + return math::equals(ca1, cb1) ? 1 + : math::equals(ca1, cb2) ? 3 + : cb1 < cb2 ? + ( ca1 < cb1 ? 0 + : ca1 > cb2 ? 4 + : 2 ) + : ( ca1 > cb1 ? 0 + : ca1 < cb2 ? 4 + : 2 ); + } }; diff --git a/boost/geometry/strategies/cartesian/centroid_average.hpp b/boost/geometry/strategies/cartesian/centroid_average.hpp index 76e2f7144c..c12f6e2024 100644 --- a/boost/geometry/strategies/cartesian/centroid_average.hpp +++ b/boost/geometry/strategies/cartesian/centroid_average.hpp @@ -4,6 +4,11 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 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 Adam Wulkiewicz, 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. @@ -15,6 +20,8 @@ #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_AVERAGE_HPP +#include <cstddef> + #include <boost/geometry/algorithms/assign.hpp> #include <boost/geometry/arithmetic/arithmetic.hpp> #include <boost/geometry/core/coordinate_type.hpp> @@ -46,7 +53,7 @@ private : class sum { friend class average; - int count; + std::size_t count; PointCentroid centroid; public : @@ -68,10 +75,15 @@ public : state.count++; } - static inline void result(sum const& state, PointCentroid& centroid) + static inline bool result(sum const& state, PointCentroid& centroid) { centroid = state.centroid; - divide_value(centroid, state.count); + if ( state.count > 0 ) + { + divide_value(centroid, state.count); + return true; + } + return false; } }; diff --git a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp index f199fb80e5..0357f17e7a 100644 --- a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp +++ b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp @@ -4,6 +4,11 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 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 Adam Wulkiewicz, 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. @@ -15,6 +20,8 @@ #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP +#include <cstddef> + #include <boost/mpl/if.hpp> #include <boost/numeric/conversion/cast.hpp> #include <boost/type_traits.hpp> @@ -143,7 +150,7 @@ private : class sums { friend class bashein_detmer; - int count; + std::size_t count; calculation_type sum_a2; calculation_type sum_x; calculation_type sum_y; 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; } diff --git a/boost/geometry/strategies/cartesian/side_of_intersection.hpp b/boost/geometry/strategies/cartesian/side_of_intersection.hpp new file mode 100644 index 0000000000..39487676c1 --- /dev/null +++ b/boost/geometry/strategies/cartesian/side_of_intersection.hpp @@ -0,0 +1,119 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP +#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP + + +#include <boost/geometry/arithmetic/determinant.hpp> +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace side +{ + +// Calculates the side of the intersection-point (if any) of +// of segment a//b w.r.t. segment c +// This is calculated without (re)calculating the IP itself again and fully +// based on integer mathematics; there are no divisions +// It can be used for either integer (rescaled) points, and also for FP +class side_of_intersection +{ +public : + + // Calculates the side of the intersection-point (if any) of + // of segment a//b w.r.t. segment c + // This is calculated without (re)calculating the IP itself again and fully + // based on integer mathematics + template <typename T, typename Segment> + static inline T side_value(Segment const& a, Segment const& b, + Segment const& c) + { + // The first point of the three segments is reused several times + T const ax = get<0, 0>(a); + T const ay = get<0, 1>(a); + T const bx = get<0, 0>(b); + T const by = get<0, 1>(b); + T const cx = get<0, 0>(c); + T const cy = get<0, 1>(c); + + T const dx_a = get<1, 0>(a) - ax; + T const dy_a = get<1, 1>(a) - ay; + + T const dx_b = get<1, 0>(b) - bx; + T const dy_b = get<1, 1>(b) - by; + + T const dx_c = get<1, 0>(c) - cx; + T const dy_c = get<1, 1>(c) - cy; + + // Cramer's rule: d (see cart_intersect.hpp) + T const d = geometry::detail::determinant<T> + ( + dx_a, dy_a, + dx_b, dy_b + ); + + T const zero = T(); + if (d == zero) + { + // There is no IP of a//b, they are collinear or parallel + // We don't have to divide but we can already conclude the side-value + // is meaningless and the resulting determinant will be 0 + return zero; + } + + // Cramer's rule: da (see cart_intersect.hpp) + T const da = geometry::detail::determinant<T> + ( + dx_b, dy_b, + ax - bx, ay - by + ); + + // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a) + // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy) + // We replace ipx by expression above and multiply each term by d + T const result = geometry::detail::determinant<T> + ( + dx_c * d, dy_c * d, + d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da + ); + + // Note: result / (d * d) + // is identical to the side_value of side_by_triangle + // Therefore, the sign is always the same as that result, and the + // resulting side (left,right,collinear) is the same + + return result; + + } + + template <typename Segment> + static inline int apply(Segment const& a, Segment const& b, Segment const& c) + { + typedef typename geometry::coordinate_type<Segment>::type coordinate_type; + coordinate_type const s = side_value<coordinate_type>(a, b, c); + coordinate_type const zero = coordinate_type(); + return math::equals(s, zero) ? 0 + : s > zero ? 1 + : -1; + } + +}; + + +}} // namespace strategy::side + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP |