diff options
Diffstat (limited to 'boost/geometry/strategies/cartesian')
10 files changed, 362 insertions, 51 deletions
diff --git a/boost/geometry/strategies/cartesian/buffer_end_round.hpp b/boost/geometry/strategies/cartesian/buffer_end_round.hpp index a233f1c4be..3d7d5bb467 100644 --- a/boost/geometry/strategies/cartesian/buffer_end_round.hpp +++ b/boost/geometry/strategies/cartesian/buffer_end_round.hpp @@ -1,6 +1,11 @@ // 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 @@ -62,8 +67,7 @@ private : DistanceType const& buffer_distance, RangeOut& range_out) const { - PromotedType const two = 2.0; - PromotedType const two_pi = two * geometry::math::pi<PromotedType>(); + PromotedType const two_pi = geometry::math::two_pi<PromotedType>(); std::size_t point_buffer_count = m_points_per_circle; diff --git a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp index 99ec80527f..5358156f6d 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp @@ -9,7 +9,7 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_MITER_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_MITER_HPP -#include <boost/assert.hpp> +#include <boost/geometry/core/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/policies/compare.hpp> #include <boost/geometry/util/math.hpp> @@ -99,7 +99,7 @@ public: if (distance > max_distance) { - BOOST_ASSERT(distance != 0.0); + BOOST_GEOMETRY_ASSERT(distance != 0.0); promoted_type const proportion = max_distance / distance; set<0>(p, get<0>(vertex) + dx * proportion); diff --git a/boost/geometry/strategies/cartesian/buffer_join_round.hpp b/boost/geometry/strategies/cartesian/buffer_join_round.hpp index 9ec51cd1ec..e1a433e1ee 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_round.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_round.hpp @@ -2,6 +2,11 @@ // 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) @@ -11,7 +16,6 @@ #include <algorithm> -#include <boost/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/policies/compare.hpp> #include <boost/geometry/strategies/buffer.hpp> @@ -19,6 +23,7 @@ #include <boost/geometry/util/select_most_precise.hpp> #ifdef BOOST_GEOMETRY_DEBUG_BUFFER_WARN +#include <iostream> #include <boost/geometry/io/wkt/wkt.hpp> #endif @@ -76,8 +81,7 @@ private : PromotedType const dx2 = get<0>(perp2) - get<0>(vertex); PromotedType const dy2 = get<1>(perp2) - get<1>(vertex); - PromotedType const two = 2.0; - PromotedType const two_pi = two * geometry::math::pi<PromotedType>(); + PromotedType const two_pi = geometry::math::two_pi<PromotedType>(); PromotedType const angle1 = atan2(dy1, dx1); PromotedType angle2 = atan2(dy2, dx2); @@ -96,14 +100,14 @@ private : // - 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 - int const n = (std::max)(static_cast<int>( - ceil(m_points_per_circle * angle_diff / two_pi)), 1); + std::size_t const n = (std::max)(static_cast<std::size_t>( + ceil(m_points_per_circle * angle_diff / two_pi)), std::size_t(1)); 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) + for (std::size_t i = 0; i < n - 1; i++, a -= diff) { Point p; set<0>(p, get<0>(vertex) + buffer_distance * cos(a)); diff --git a/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp b/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp index 1444c795af..ed3a010cd5 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp @@ -9,7 +9,6 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_BY_DIVIDE_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_BY_DIVIDE_HPP -#include <boost/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/policies/compare.hpp> #include <boost/geometry/strategies/buffer.hpp> diff --git a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp index 86ebc43c9c..f289857177 100644 --- a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp +++ b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp @@ -85,8 +85,7 @@ public : promoted_type const buffer_distance = distance_strategy.apply(point, point, strategy::buffer::buffer_side_left); - promoted_type const two = 2.0; - promoted_type const two_pi = two * geometry::math::pi<promoted_type>(); + promoted_type const two_pi = geometry::math::two_pi<promoted_type>(); promoted_type const diff = two_pi / promoted_type(m_count); promoted_type a = 0; diff --git a/boost/geometry/strategies/cartesian/buffer_side_straight.hpp b/boost/geometry/strategies/cartesian/buffer_side_straight.hpp index 24655ab3d7..75247377fa 100644 --- a/boost/geometry/strategies/cartesian/buffer_side_straight.hpp +++ b/boost/geometry/strategies/cartesian/buffer_side_straight.hpp @@ -9,6 +9,8 @@ #include <cstddef> +#include <boost/math/special_functions/fpclassify.hpp> + #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/access.hpp> #include <boost/geometry/util/math.hpp> @@ -51,9 +53,9 @@ public : typename OutputRange, typename DistanceStrategy > - static inline void apply( + static inline result_code apply( Point const& input_p1, Point const& input_p2, - strategy::buffer::buffer_side_selector side, + buffer_side_selector side, DistanceStrategy const& distance, OutputRange& output_range) { @@ -73,19 +75,46 @@ public : // For normalization [0,1] (=dot product d.d, sqrt) promoted_type const length = geometry::math::sqrt(dx * dx + dy * dy); + if (! boost::math::isfinite(length)) + { + // In case of coordinates differences of e.g. 1e300, length + // will overflow and we should not generate output +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_WARN + std::cout << "Error in length calculation for points " + << geometry::wkt(input_p1) << " " << geometry::wkt(input_p2) + << " length: " << length << std::endl; +#endif + return result_error_numerical; + } + if (geometry::math::equals(length, 0)) { // Coordinates are simplified and therefore most often not equal. // But if simplify is skipped, or for lines with two // equal points, length is 0 and we cannot generate output. - return; + return result_no_output; } + promoted_type const d = distance.apply(input_p1, input_p2, side); + // Generate the normalized perpendicular p, to the left (ccw) promoted_type const px = -dy / length; promoted_type const py = dx / length; - promoted_type const d = distance.apply(input_p1, input_p2, side); + if (geometry::math::equals(px, 0) + && geometry::math::equals(py, 0)) + { + // This basically should not occur - because of the checks above. + // There are no unit tests triggering this condition +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_WARN + std::cout << "Error in perpendicular calculation for points " + << geometry::wkt(input_p1) << " " << geometry::wkt(input_p2) + << " length: " << length + << " distance: " << d + << std::endl; +#endif + return result_no_output; + } output_range.resize(2); @@ -93,6 +122,8 @@ public : set<1>(output_range.front(), get<1>(input_p1) + py * d); set<0>(output_range.back(), get<0>(input_p2) + px * d); set<1>(output_range.back(), get<1>(input_p2) + py * d); + + return result_normal; } #endif // DOXYGEN_SHOULD_SKIP_THIS }; diff --git a/boost/geometry/strategies/cartesian/cart_intersect.hpp b/boost/geometry/strategies/cartesian/cart_intersect.hpp index a7bd385226..8a9857376a 100644 --- a/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -119,15 +119,10 @@ struct relate_cartesian_segments boost::ignore_unused_variable_warning(robust_policy); - typedef typename select_calculation_type - <Segment1, Segment2, CalculationType>::type coordinate_type; - using geometry::detail::equals::equals_point_point; bool const a_is_point = equals_point_point(robust_a1, robust_a2); bool const b_is_point = equals_point_point(robust_b1, robust_b2); - typedef side::side_by_triangle<coordinate_type> side; - if(a_is_point && b_is_point) { return equals_point_point(robust_a1, robust_b2) @@ -136,20 +131,32 @@ struct relate_cartesian_segments ; } + typedef typename select_calculation_type + <Segment1, Segment2, CalculationType>::type coordinate_type; + + typedef side::side_by_triangle<coordinate_type> side; + side_info sides; sides.set<0>(side::apply(robust_b1, robust_b2, robust_a1), side::apply(robust_b1, robust_b2, robust_a2)); - sides.set<1>(side::apply(robust_a1, robust_a2, robust_b1), - side::apply(robust_a1, robust_a2, robust_b2)); - bool collinear = sides.collinear(); + if (sides.same<0>()) + { + // Both points are at same side of other segment, we can leave + return Policy::disjoint(); + } - if (sides.same<0>() || sides.same<1>()) + sides.set<1>(side::apply(robust_a1, robust_a2, robust_b1), + side::apply(robust_a1, robust_a2, robust_b2)); + + if (sides.same<1>()) { // Both points are at same side of other segment, we can leave return Policy::disjoint(); } + bool collinear = sides.collinear(); + typedef typename select_most_precise < coordinate_type, double @@ -264,7 +271,7 @@ private: 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"); + //BOOST_GEOMETRY_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 diff --git a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp index 0357f17e7a..47763a9f30 100644 --- a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp +++ b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp @@ -1,8 +1,8 @@ // 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. @@ -22,6 +22,7 @@ #include <cstddef> +#include <boost/math/special_functions/fpclassify.hpp> #include <boost/mpl/if.hpp> #include <boost/numeric/conversion/cast.hpp> #include <boost/type_traits.hpp> @@ -206,11 +207,18 @@ public : Point >::type coordinate_type; - set<0>(centroid, - boost::numeric_cast<coordinate_type>(state.sum_x / a3)); - set<1>(centroid, - boost::numeric_cast<coordinate_type>(state.sum_y / a3)); - return true; + // Prevent NaN centroid coordinates + if (boost::math::isfinite(a3)) + { + // NOTE: above calculation_type is checked, not the centroid coordinate_type + // which means that the centroid can still be filled with INF + // if e.g. calculation_type is double and centroid contains floats + set<0>(centroid, + boost::numeric_cast<coordinate_type>(state.sum_x / a3)); + set<1>(centroid, + boost::numeric_cast<coordinate_type>(state.sum_y / a3)); + return true; + } } return false; diff --git a/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp b/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp index b788738d15..0735e925b6 100644 --- a/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp +++ b/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp @@ -1,7 +1,12 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2009-2012 Mateusz Loskot, London, UK. -// Copyright (c) 2009-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. +// Copyright (c) 2009-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 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. @@ -13,9 +18,13 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_WEIGHTED_LENGTH_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_WEIGHTED_LENGTH_HPP +#include <boost/math/special_functions/fpclassify.hpp> +#include <boost/numeric/conversion/cast.hpp> + #include <boost/geometry/algorithms/detail/distance/interface.hpp> #include <boost/geometry/algorithms/detail/distance/point_to_geometry.hpp> #include <boost/geometry/arithmetic/arithmetic.hpp> +#include <boost/geometry/util/for_each_coordinate.hpp> #include <boost/geometry/util/select_most_precise.hpp> #include <boost/geometry/strategies/centroid.hpp> #include <boost/geometry/strategies/default_distance_result.hpp> @@ -91,17 +100,37 @@ public : static inline bool result(state_type const& state, Point& centroid) { distance_type const zero = distance_type(); - if (! geometry::math::equals(state.length, zero)) + if (! geometry::math::equals(state.length, zero) + && boost::math::isfinite(state.length)) // Prevent NaN centroid coordinates { - assign_zero(centroid); - add_point(centroid, state.average_sum); - divide_value(centroid, state.length); + // NOTE: above distance_type is checked, not the centroid coordinate_type + // which means that the centroid can still be filled with INF + // if e.g. distance_type is double and centroid contains floats + geometry::for_each_coordinate(centroid, set_sum_div_length(state)); return true; } return false; } + struct set_sum_div_length + { + state_type const& m_state; + set_sum_div_length(state_type const& state) + : m_state(state) + {} + template <typename Pt, std::size_t Dimension> + void apply(Pt & centroid) const + { + typedef typename geometry::coordinate_type<Pt>::type coordinate_type; + geometry::set<Dimension>( + centroid, + boost::numeric_cast<coordinate_type>( + geometry::get<Dimension>(m_state.average_sum) / m_state.length + ) + ); + } + }; }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS diff --git a/boost/geometry/strategies/cartesian/side_of_intersection.hpp b/boost/geometry/strategies/cartesian/side_of_intersection.hpp index 39487676c1..db57644ad5 100644 --- a/boost/geometry/strategies/cartesian/side_of_intersection.hpp +++ b/boost/geometry/strategies/cartesian/side_of_intersection.hpp @@ -2,6 +2,11 @@ // Copyright (c) 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 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 // http://www.boost.org/LICENSE_1_0.txt) @@ -10,11 +15,25 @@ #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP +#include <limits> + +#include <boost/core/ignore_unused.hpp> +#include <boost/type_traits/is_integral.hpp> +#include <boost/type_traits/make_unsigned.hpp> + #include <boost/geometry/arithmetic/determinant.hpp> #include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/assert.hpp> #include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> +#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> #include <boost/geometry/util/math.hpp> +#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG +#include <boost/math/common_factor_ct.hpp> +#include <boost/math/common_factor_rt.hpp> +#include <boost/multiprecision/cpp_int.hpp> +#endif namespace boost { namespace geometry { @@ -22,6 +41,94 @@ namespace boost { namespace geometry namespace strategy { namespace side { +namespace detail +{ + +// A tool for multiplication of integers avoiding overflow +// It's a temporary workaround until we can use Multiprecision +// The algorithm is based on Karatsuba algorithm +// see: http://en.wikipedia.org/wiki/Karatsuba_algorithm +template <typename T> +struct multiplicable_integral +{ + // Currently this tool can't be used with non-integral coordinate types. + // Also side_of_intersection strategy sign_of_product() and sign_of_compare() + // functions would have to be modified to properly support floating-point + // types (comparisons and multiplication). + BOOST_STATIC_ASSERT(boost::is_integral<T>::value); + + static const std::size_t bits = CHAR_BIT * sizeof(T); + static const std::size_t half_bits = bits / 2; + typedef typename boost::make_unsigned<T>::type unsigned_type; + static const unsigned_type base = unsigned_type(1) << half_bits; // 2^half_bits + + int m_sign; + unsigned_type m_ms; + unsigned_type m_ls; + + multiplicable_integral(int sign, unsigned_type ms, unsigned_type ls) + : m_sign(sign), m_ms(ms), m_ls(ls) + {} + + explicit multiplicable_integral(T const& val) + { + unsigned_type val_u = val > 0 ? + unsigned_type(val) + : val == (std::numeric_limits<T>::min)() ? + unsigned_type((std::numeric_limits<T>::max)()) + 1 + : unsigned_type(-val); + // MMLL -> S 00MM 00LL + m_sign = math::sign(val); + m_ms = val_u >> half_bits; // val_u / base + m_ls = val_u - m_ms * base; + } + + friend multiplicable_integral operator*(multiplicable_integral const& a, + multiplicable_integral const& b) + { + // (S 00MM 00LL) * (S 00MM 00LL) -> (S Z2MM 00LL) + unsigned_type z2 = a.m_ms * b.m_ms; + unsigned_type z0 = a.m_ls * b.m_ls; + unsigned_type z1 = (a.m_ms + a.m_ls) * (b.m_ms + b.m_ls) - z2 - z0; + // z0 may be >= base so it must be normalized to allow comparison + unsigned_type z0_ms = z0 >> half_bits; // z0 / base + return multiplicable_integral(a.m_sign * b.m_sign, + z2 * base + z1 + z0_ms, + z0 - base * z0_ms); + } + + friend bool operator<(multiplicable_integral const& a, + multiplicable_integral const& b) + { + if ( a.m_sign == b.m_sign ) + { + bool u_less = a.m_ms < b.m_ms + || (a.m_ms == b.m_ms && a.m_ls < b.m_ls); + return a.m_sign > 0 ? u_less : (! u_less); + } + else + { + return a.m_sign < b.m_sign; + } + } + + friend bool operator>(multiplicable_integral const& a, + multiplicable_integral const& b) + { + return b < a; + } + + template <typename CmpVal> + void check_value(CmpVal const& cmp_val) const + { + unsigned_type b = base; // a workaround for MinGW - undefined reference base + CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls)); + BOOST_GEOMETRY_ASSERT(cmp_val == val); + } +}; + +} // namespace detail + // 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 @@ -29,15 +136,93 @@ namespace strategy { namespace side // It can be used for either integer (rescaled) points, and also for FP class side_of_intersection { +private : + template <typename T, typename U> + static inline + int sign_of_product(T const& a, U const& b) + { + return a == 0 || b == 0 ? 0 + : a > 0 && b > 0 ? 1 + : a < 0 && b < 0 ? 1 + : -1; + } + + template <typename T> + static inline + int sign_of_compare(T const& a, T const& b, T const& c, T const& d) + { + // Both a*b and c*d are positive + // We have to judge if a*b > c*d + + using side::detail::multiplicable_integral; + multiplicable_integral<T> ab = multiplicable_integral<T>(a) + * multiplicable_integral<T>(b); + multiplicable_integral<T> cd = multiplicable_integral<T>(c) + * multiplicable_integral<T>(d); + + int result = ab > cd ? 1 + : ab < cd ? -1 + : 0 + ; + +#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG + using namespace boost::multiprecision; + cpp_int const lab = cpp_int(a) * cpp_int(b); + cpp_int const lcd = cpp_int(c) * cpp_int(d); + + ab.check_value(lab); + cd.check_value(lcd); + + int result2 = lab > lcd ? 1 + : lab < lcd ? -1 + : 0 + ; + BOOST_GEOMETRY_ASSERT(result == result2); +#endif + + return result; + } + + template <typename T> + static inline + int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d) + { + // sign of a*b+c*d, 1 if positive, -1 if negative, else 0 + int const ab = sign_of_product(a, b); + int const cd = sign_of_product(c, d); + if (ab == 0) + { + return cd; + } + if (cd == 0) + { + return ab; + } + + if (ab == cd) + { + // Both positive or both negative + return ab; + } + + // One is positive, one is negative, both are non zero + // If ab is positive, we have to judge if a*b > -c*d (then 1 because sum is positive) + // If ab is negative, we have to judge if c*d > -a*b (idem) + return ab == 1 + ? sign_of_compare(a, b, -c, d) + : sign_of_compare(c, d, -a, b); + } + + 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> + template <typename T, typename Segment, typename Point> static inline T side_value(Segment const& a, Segment const& b, - Segment const& c) + Segment const& c, Point const& fallback_point) { // The first point of the three segments is reused several times T const ax = get<0, 0>(a); @@ -67,9 +252,15 @@ public : 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; + // Assuming they intersect (this method should be called for + // segments known to intersect), they are collinear and overlap. + // They have one or two intersection points - we don't know and + // have to rely on the fallback intersection point + + Point c1, c2; + geometry::detail::assign_point_from_index<0>(c, c1); + geometry::detail::assign_point_from_index<1>(c, c2); + return side_by_triangle<>::apply(c1, c2, fallback_point); } // Cramer's rule: da (see cart_intersect.hpp) @@ -82,7 +273,9 @@ public : // 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> + +#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG + T const result1 = geometry::detail::determinant<T> ( dx_c * d, dy_c * d, d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da @@ -93,15 +286,52 @@ public : // Therefore, the sign is always the same as that result, and the // resulting side (left,right,collinear) is the same + // The first row we divide again by d because of determinant multiply rule + T const result2 = d * geometry::detail::determinant<T> + ( + dx_c, dy_c, + d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da + ); + // Write out: + T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da) + - dy_c * (d * (ax - cx) + dx_a * da)); + // Write out in braces: + T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da + - dy_c * d * (ax - cx) - dy_c * dx_a * da); + // Write in terms of d * XX + da * YY + T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx)) + + da * (dx_c * dy_a - dy_c * dx_a)); + + boost::ignore_unused(result1, result2, result3, result4, result5); + //return result; +#endif + + // We consider the results separately + // (in the end we only have to return the side-value 1,0 or -1) + + // To avoid multiplications we judge the product (easy, avoids *d) + // and the sign of p*q+r*s (more elaborate) + T const result = sign_of_product + ( + d, + sign_of_addition_of_two_products + ( + d, dx_c * (ay - cy) - dy_c * (ax - cx), + da, dx_c * dy_a - dy_c * dx_a + ) + ); return result; + } - template <typename Segment> - static inline int apply(Segment const& a, Segment const& b, Segment const& c) + template <typename Segment, typename Point> + static inline int apply(Segment const& a, Segment const& b, + Segment const& c, + Point const& fallback_point) { typedef typename geometry::coordinate_type<Segment>::type coordinate_type; - coordinate_type const s = side_value<coordinate_type>(a, b, c); + coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point); coordinate_type const zero = coordinate_type(); return math::equals(s, zero) ? 0 : s > zero ? 1 |