Chris@102: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@102: 
Chris@102: // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
Chris@102: 
Chris@102: // Use, modification and distribution is subject to the Boost Software License,
Chris@102: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@102: // http://www.boost.org/LICENSE_1_0.txt)
Chris@102: 
Chris@102: #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
Chris@102: #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
Chris@102: 
Chris@102: 
Chris@102: #include <boost/geometry/arithmetic/determinant.hpp>
Chris@102: #include <boost/geometry/core/access.hpp>
Chris@102: #include <boost/geometry/core/coordinate_type.hpp>
Chris@102: #include <boost/geometry/util/math.hpp>
Chris@102: 
Chris@102: 
Chris@102: namespace boost { namespace geometry
Chris@102: {
Chris@102: 
Chris@102: namespace strategy { namespace side
Chris@102: {
Chris@102: 
Chris@102: // Calculates the side of the intersection-point (if any) of
Chris@102: // of segment a//b w.r.t. segment c
Chris@102: // This is calculated without (re)calculating the IP itself again and fully
Chris@102: // based on integer mathematics; there are no divisions
Chris@102: // It can be used for either integer (rescaled) points, and also for FP
Chris@102: class side_of_intersection
Chris@102: {
Chris@102: public :
Chris@102: 
Chris@102:     // Calculates the side of the intersection-point (if any) of
Chris@102:     // of segment a//b w.r.t. segment c
Chris@102:     // This is calculated without (re)calculating the IP itself again and fully
Chris@102:     // based on integer mathematics
Chris@102:     template <typename T, typename Segment>
Chris@102:     static inline T side_value(Segment const& a, Segment const& b,
Chris@102:                 Segment const& c)
Chris@102:     {
Chris@102:         // The first point of the three segments is reused several times
Chris@102:         T const ax = get<0, 0>(a);
Chris@102:         T const ay = get<0, 1>(a);
Chris@102:         T const bx = get<0, 0>(b);
Chris@102:         T const by = get<0, 1>(b);
Chris@102:         T const cx = get<0, 0>(c);
Chris@102:         T const cy = get<0, 1>(c);
Chris@102: 
Chris@102:         T const dx_a = get<1, 0>(a) - ax;
Chris@102:         T const dy_a = get<1, 1>(a) - ay;
Chris@102: 
Chris@102:         T const dx_b = get<1, 0>(b) - bx;
Chris@102:         T const dy_b = get<1, 1>(b) - by;
Chris@102: 
Chris@102:         T const dx_c = get<1, 0>(c) - cx;
Chris@102:         T const dy_c = get<1, 1>(c) - cy;
Chris@102: 
Chris@102:         // Cramer's rule: d (see cart_intersect.hpp)
Chris@102:         T const d = geometry::detail::determinant<T>
Chris@102:                     (
Chris@102:                         dx_a, dy_a,
Chris@102:                         dx_b, dy_b
Chris@102:                     );
Chris@102: 
Chris@102:         T const zero = T();
Chris@102:         if (d == zero)
Chris@102:         {
Chris@102:             // There is no IP of a//b, they are collinear or parallel
Chris@102:             // We don't have to divide but we can already conclude the side-value
Chris@102:             // is meaningless and the resulting determinant will be 0
Chris@102:             return zero;
Chris@102:         }
Chris@102: 
Chris@102:         // Cramer's rule: da (see cart_intersect.hpp)
Chris@102:         T const da = geometry::detail::determinant<T>
Chris@102:                     (
Chris@102:                         dx_b,    dy_b,
Chris@102:                         ax - bx, ay - by
Chris@102:                     );
Chris@102: 
Chris@102:         // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
Chris@102:         // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
Chris@102:         // We replace ipx by expression above and multiply each term by d
Chris@102:         T const result = geometry::detail::determinant<T>
Chris@102:                     (
Chris@102:                         dx_c * d,                   dy_c * d,
Chris@102:                         d * (ax - cx) + dx_a * da,  d * (ay - cy) + dy_a * da
Chris@102:                     );
Chris@102: 
Chris@102:         // Note: result / (d * d)
Chris@102:         // is identical to the side_value of side_by_triangle
Chris@102:         // Therefore, the sign is always the same as that result, and the
Chris@102:         // resulting side (left,right,collinear) is the same
Chris@102: 
Chris@102:         return result;
Chris@102: 
Chris@102:     }
Chris@102: 
Chris@102:     template <typename Segment>
Chris@102:     static inline int apply(Segment const& a, Segment const& b, Segment const& c)
Chris@102:     {
Chris@102:         typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
Chris@102:         coordinate_type const s = side_value<coordinate_type>(a, b, c);
Chris@102:         coordinate_type const zero = coordinate_type();
Chris@102:         return math::equals(s, zero) ? 0
Chris@102:             : s > zero ? 1
Chris@102:             : -1;
Chris@102:     }
Chris@102: 
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: }} // namespace strategy::side
Chris@102: 
Chris@102: }} // namespace boost::geometry
Chris@102: 
Chris@102: 
Chris@102: #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP