Chris@16: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@16: 
Chris@16: // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
Chris@16: 
Chris@16: // Use, modification and distribution is subject to the Boost Software License,
Chris@16: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_GEOMETRY_GEOMETRY_POLICIES_RELATE_INTERSECTION_POINTS_HPP
Chris@16: #define BOOST_GEOMETRY_GEOMETRY_POLICIES_RELATE_INTERSECTION_POINTS_HPP
Chris@16: 
Chris@16: 
Chris@16: #include <algorithm>
Chris@16: #include <string>
Chris@16: 
Chris@16: #include <boost/concept_check.hpp>
Chris@16: #include <boost/numeric/conversion/cast.hpp>
Chris@16: 
Chris@101: #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
Chris@16: #include <boost/geometry/core/access.hpp>
Chris@16: #include <boost/geometry/strategies/side_info.hpp>
Chris@101: #include <boost/geometry/util/promote_integral.hpp>
Chris@16: #include <boost/geometry/util/select_calculation_type.hpp>
Chris@16: #include <boost/geometry/util/select_most_precise.hpp>
Chris@101: #include <boost/geometry/util/math.hpp>
Chris@16: 
Chris@16: namespace boost { namespace geometry
Chris@16: {
Chris@16: 
Chris@16: namespace policies { namespace relate
Chris@16: {
Chris@16: 
Chris@16: 
Chris@101: /*!
Chris@101: \brief Policy calculating the intersection points themselves
Chris@101:  */
Chris@101: template
Chris@101: <
Chris@101:     typename ReturnType
Chris@101: >
Chris@16: struct segments_intersection_points
Chris@16: {
Chris@16:     typedef ReturnType return_type;
Chris@16: 
Chris@101:     template
Chris@101:     <
Chris@101:         typename Point,
Chris@101:         typename Segment,
Chris@101:         typename SegmentRatio,
Chris@101:         typename T
Chris@101:     >
Chris@101:     static inline void assign(Point& point,
Chris@101:                 Segment const& segment,
Chris@101:                 SegmentRatio const& ratio,
Chris@101:                 T const& dx, T const& dy)
Chris@101:     {
Chris@101:         typedef typename geometry::coordinate_type<Point>::type coordinate_type;
Chris@16: 
Chris@101:         // Calculate the intersection point based on segment_ratio
Chris@101:         // Up to now, division was postponed. Here we divide using numerator/
Chris@101:         // denominator. In case of integer this results in an integer
Chris@101:         // division.
Chris@101:         BOOST_ASSERT(ratio.denominator() != 0);
Chris@101: 
Chris@101:         typedef typename promote_integral<coordinate_type>::type promoted_type;
Chris@101: 
Chris@101:         promoted_type const numerator
Chris@101:             = boost::numeric_cast<promoted_type>(ratio.numerator());
Chris@101:         promoted_type const denominator
Chris@101:             = boost::numeric_cast<promoted_type>(ratio.denominator());
Chris@101:         promoted_type const dx_promoted = boost::numeric_cast<promoted_type>(dx);
Chris@101:         promoted_type const dy_promoted = boost::numeric_cast<promoted_type>(dy);
Chris@101: 
Chris@101:         set<0>(point, get<0, 0>(segment) + boost::numeric_cast
Chris@101:             <
Chris@101:                 coordinate_type
Chris@101:             >(numerator * dx_promoted / denominator));
Chris@101:         set<1>(point, get<0, 1>(segment) + boost::numeric_cast
Chris@101:             <
Chris@101:                 coordinate_type
Chris@101:             >(numerator * dy_promoted / denominator));
Chris@101:     }
Chris@101: 
Chris@101: 
Chris@101:     template
Chris@101:     <
Chris@101:         typename Segment1,
Chris@101:         typename Segment2,
Chris@101:         typename SegmentIntersectionInfo
Chris@101:     >
Chris@101:     static inline return_type segments_crosses(side_info const&,
Chris@101:                     SegmentIntersectionInfo const& sinfo,
Chris@101:                     Segment1 const& s1, Segment2 const& s2)
Chris@16:     {
Chris@16:         return_type result;
Chris@16:         result.count = 1;
Chris@101: 
Chris@101:         if (sinfo.robust_ra < sinfo.robust_rb)
Chris@101:         {
Chris@101:             assign(result.intersections[0], s1, sinfo.robust_ra,
Chris@101:                 sinfo.dx_a, sinfo.dy_a);
Chris@101:         }
Chris@101:         else
Chris@101:         {
Chris@101:             assign(result.intersections[0], s2, sinfo.robust_rb,
Chris@101:                 sinfo.dx_b, sinfo.dy_b);
Chris@101:         }
Chris@101: 
Chris@101:         result.fractions[0].assign(sinfo);
Chris@16: 
Chris@16:         return result;
Chris@16:     }
Chris@16: 
Chris@101:     template <typename Segment1, typename Segment2, typename Ratio>
Chris@101:     static inline return_type segments_collinear(
Chris@101:         Segment1 const& a, Segment2 const& b, bool /*opposite*/,
Chris@101:         int a1_wrt_b, int a2_wrt_b, int b1_wrt_a, int b2_wrt_a,
Chris@101:         Ratio const& ra_from_wrt_b, Ratio const& ra_to_wrt_b,
Chris@101:         Ratio const& rb_from_wrt_a, Ratio const& rb_to_wrt_a)
Chris@16:     {
Chris@16:         return_type result;
Chris@101:         unsigned int index = 0, count_a = 0, count_b = 0;
Chris@101:         Ratio on_a[2];
Chris@16: 
Chris@101:         // The conditions "index < 2" are necessary for non-robust handling,
Chris@101:         // if index would be 2 this indicate an (currently uncatched) error
Chris@16: 
Chris@101:         // IMPORTANT: the order of conditions is different as in direction.hpp
Chris@101:         if (a1_wrt_b >= 1 && a1_wrt_b <= 3 // ra_from_wrt_b.on_segment()
Chris@101:             && index < 2)
Chris@101:         {
Chris@101:             //     a1--------->a2
Chris@101:             // b1----->b2
Chris@101:             //
Chris@101:             // ra1 (relative to b) is between 0/1:
Chris@101:             // -> First point of A is intersection point
Chris@101:             detail::assign_point_from_index<0>(a, result.intersections[index]);
Chris@101:             result.fractions[index].assign(Ratio::zero(), ra_from_wrt_b);
Chris@101:             on_a[index] = Ratio::zero();
Chris@101:             index++;
Chris@101:             count_a++;
Chris@101:         }
Chris@101:         if (b1_wrt_a == 2 //rb_from_wrt_a.in_segment()
Chris@101:             && index < 2)
Chris@101:         {
Chris@101:             // We take the first intersection point of B
Chris@101:             // a1--------->a2
Chris@101:             //         b1----->b2
Chris@101:             // But only if it is not located on A
Chris@101:             // a1--------->a2
Chris@101:             // b1----->b2      rb_from_wrt_a == 0/1 -> a already taken
Chris@16: 
Chris@101:             detail::assign_point_from_index<0>(b, result.intersections[index]);
Chris@101:             result.fractions[index].assign(rb_from_wrt_a, Ratio::zero());
Chris@101:             on_a[index] = rb_from_wrt_a;
Chris@101:             index++;
Chris@101:             count_b++;
Chris@101:         }
Chris@16: 
Chris@101:         if (a2_wrt_b >= 1 && a2_wrt_b <= 3 //ra_to_wrt_b.on_segment()
Chris@101:             && index < 2)
Chris@101:         {
Chris@101:             // Similarly, second IP (here a2)
Chris@101:             // a1--------->a2
Chris@101:             //         b1----->b2
Chris@101:             detail::assign_point_from_index<1>(a, result.intersections[index]);
Chris@101:             result.fractions[index].assign(Ratio::one(), ra_to_wrt_b);
Chris@101:             on_a[index] = Ratio::one();
Chris@101:             index++;
Chris@101:             count_a++;
Chris@101:         }
Chris@101:         if (b2_wrt_a == 2 // rb_to_wrt_a.in_segment()
Chris@101:             && index < 2)
Chris@101:         {
Chris@101:             detail::assign_point_from_index<1>(b, result.intersections[index]);
Chris@101:             result.fractions[index].assign(rb_to_wrt_a, Ratio::one());
Chris@101:             on_a[index] = rb_to_wrt_a;
Chris@101:             index++;
Chris@101:             count_b++;
Chris@101:         }
Chris@16: 
Chris@101:         // TEMPORARY
Chris@101:         // If both are from b, and b is reversed w.r.t. a, we swap IP's
Chris@101:         // to align them w.r.t. a
Chris@101:         // get_turn_info still relies on some order (in some collinear cases)
Chris@101:         if (index == 2 && on_a[1] < on_a[0])
Chris@101:         {
Chris@101:             std::swap(result.fractions[0], result.fractions[1]);
Chris@101:             std::swap(result.intersections[0], result.intersections[1]);
Chris@101:         }
Chris@101: 
Chris@101:         result.count = index;
Chris@101: 
Chris@16:         return result;
Chris@16:     }
Chris@16: 
Chris@16:     static inline return_type disjoint()
Chris@16:     {
Chris@16:         return return_type();
Chris@16:     }
Chris@16:     static inline return_type error(std::string const&)
Chris@16:     {
Chris@16:         return return_type();
Chris@16:     }
Chris@16: 
Chris@101:     // Both degenerate
Chris@101:     template <typename Segment>
Chris@101:     static inline return_type degenerate(Segment const& segment, bool)
Chris@16:     {
Chris@16:         return_type result;
Chris@16:         result.count = 1;
Chris@101:         set<0>(result.intersections[0], get<0, 0>(segment));
Chris@101:         set<1>(result.intersections[0], get<0, 1>(segment));
Chris@101:         return result;
Chris@101:     }
Chris@101: 
Chris@101:     // One degenerate
Chris@101:     template <typename Segment, typename Ratio>
Chris@101:     static inline return_type one_degenerate(Segment const& degenerate_segment,
Chris@101:             Ratio const& ratio, bool a_degenerate)
Chris@101:     {
Chris@101:         return_type result;
Chris@101:         result.count = 1;
Chris@101:         set<0>(result.intersections[0], get<0, 0>(degenerate_segment));
Chris@101:         set<1>(result.intersections[0], get<0, 1>(degenerate_segment));
Chris@101:         if (a_degenerate)
Chris@101:         {
Chris@101:             // IP lies on ratio w.r.t. segment b
Chris@101:             result.fractions[0].assign(Ratio::zero(), ratio);
Chris@101:         }
Chris@101:         else
Chris@101:         {
Chris@101:             result.fractions[0].assign(ratio, Ratio::zero());
Chris@101:         }
Chris@16:         return result;
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: }} // namespace policies::relate
Chris@16: 
Chris@16: }} // namespace boost::geometry
Chris@16: 
Chris@16: #endif // BOOST_GEOMETRY_GEOMETRY_POLICIES_RELATE_INTERSECTION_POINTS_HPP