detail/buffer/line_line_intersection.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2012-2020 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // Use, modification and distribution is subject to the Boost Software License,
0006 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0007 // http://www.boost.org/LICENSE_1_0.txt)
0008
0009 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_LINE_LINE_INTERSECTION_HPP
0010 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_LINE_LINE_INTERSECTION_HPP
0011
0012 #include <boost/geometry/algorithms/detail/make/make.hpp>
0013 #include <boost/geometry/arithmetic/infinite_line_functions.hpp>
0014 #include <boost/geometry/util/math.hpp>
0015
0016 namespace boost { namespace geometry
0017 {
0018
0019
0020 #ifndef DOXYGEN_NO_DETAIL
0021 namespace detail { namespace buffer
0022 {
0023
0024 struct line_line_intersection
0025 {
0026     template <typename Point>
0027     static Point between_point(Point const& a, Point const& b)
0028     {
0029         Point result;
0030         geometry::set<0>(result, (geometry::get<0>(a) + geometry::get<0>(b)) / 2.0);
0031         geometry::set<1>(result, (geometry::get<1>(a) + geometry::get<1>(b)) / 2.0);
0032         return result;
0033     }
0034
0035     template <typename Point>
0036     static bool
0037     apply(Point const& pi, Point const& pj, Point const& qi, Point const& qj,
0038           Point const& vertex, bool equidistant, Point& ip)
0039     {
0040         // Calculates ip (below) by either intersecting p (pi, pj)
0041         // with q (qi, qj) or by taking a point between pj and qi (b) and
0042         // intersecting r (b, v), where v is the original vertex, with p (or q).
0043         // The reason for dual approach: p might be nearly collinear with q,
0044         // and in that case the intersection points can lose precision
0045         // (or be plainly wrong).
0046         // Therefore it takes the most precise option (this is usually p, r)
0047         //
0048         //             /qj                     |
0049         //            /                        |
0050         //           /      /                  |
0051         //          /      /                   |
0052         //         /      /                    |
0053         //        /qi    /                     |
0054         //              /                      |
0055         //   ip *  + b * v                     |
0056         //              \                      |
0057         //        \pj    \                     |
0058         //         \      \                    |
0059         //          \      \                   |
0060         //           \      \                  |
0061         //            \pi    \                 |
0062         //
0063         // If generated sides along the segments can have an adapted distance,
0064         // in a custom strategy, then the calculation of the point in between
0065         // might be incorrect and the optimization is not used.
0066
0067         using ct = typename coordinate_type<Point>::type;
0068
0069         auto const p = detail::make::make_infinite_line<ct>(pi, pj);
0070         auto const q = detail::make::make_infinite_line<ct>(qi, qj);
0071
0072         using line = decltype(p);
0073         using arithmetic::determinant;
0074         using arithmetic::assign_intersection_point;
0075
0076         // The denominator is the determinant of (a,b) values of lines p q
0077         // | pa pa |
0078         // | qb qb |
0079         auto const denominator_pq = determinant<line, &line::a, &line::b>(p, q);
0080         static decltype(denominator_pq) const zero = 0;
0081
0082         if (equidistant)
0083         {
0084             auto const between = between_point(pj, qi);
0085             auto const r = detail::make::make_infinite_line<ct>(vertex, between);
0086             auto const denominator_pr = determinant<line, &line::a, &line::b>(p, r);
0087
0088             if (math::equals(denominator_pq, zero)
0089                 && math::equals(denominator_pr, zero))
0090             {
0091                 // Degenerate case (for example when length results in <inf>)
0092                 return false;
0093             }
0094
0095             ip = geometry::math::abs(denominator_pq) > geometry::math::abs(denominator_pr)
0096                  ? assign_intersection_point<Point>(p, q, denominator_pq)
0097                  : assign_intersection_point<Point>(p, r, denominator_pr);
0098         }
0099         else
0100         {
0101             if (math::equals(denominator_pq, zero))
0102             {
0103                 return false;
0104             }
0105             ip = assign_intersection_point<Point>(p, q, denominator_pq);
0106         }
0107
0108         return true;
0109     }
0110 };
0111
0112
0113 }} // namespace detail::buffer
0114 #endif // DOXYGEN_NO_DETAIL
0115
0116
0117 }} // namespace boost::geometry
0118
0119
0120 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_LINE_LINE_INTERSECTION_HPP