Chris@102
|
1 // Boost.Geometry (aka GGL, Generic Geometry Library)
|
Chris@102
|
2
|
Chris@102
|
3 // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
|
Chris@102
|
4
|
Chris@102
|
5 // Use, modification and distribution is subject to the Boost Software License,
|
Chris@102
|
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
Chris@102
|
7 // http://www.boost.org/LICENSE_1_0.txt)
|
Chris@102
|
8
|
Chris@102
|
9 #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
|
Chris@102
|
10 #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
|
Chris@102
|
11
|
Chris@102
|
12
|
Chris@102
|
13 #include <boost/geometry/arithmetic/determinant.hpp>
|
Chris@102
|
14 #include <boost/geometry/core/access.hpp>
|
Chris@102
|
15 #include <boost/geometry/core/coordinate_type.hpp>
|
Chris@102
|
16 #include <boost/geometry/util/math.hpp>
|
Chris@102
|
17
|
Chris@102
|
18
|
Chris@102
|
19 namespace boost { namespace geometry
|
Chris@102
|
20 {
|
Chris@102
|
21
|
Chris@102
|
22 namespace strategy { namespace side
|
Chris@102
|
23 {
|
Chris@102
|
24
|
Chris@102
|
25 // Calculates the side of the intersection-point (if any) of
|
Chris@102
|
26 // of segment a//b w.r.t. segment c
|
Chris@102
|
27 // This is calculated without (re)calculating the IP itself again and fully
|
Chris@102
|
28 // based on integer mathematics; there are no divisions
|
Chris@102
|
29 // It can be used for either integer (rescaled) points, and also for FP
|
Chris@102
|
30 class side_of_intersection
|
Chris@102
|
31 {
|
Chris@102
|
32 public :
|
Chris@102
|
33
|
Chris@102
|
34 // Calculates the side of the intersection-point (if any) of
|
Chris@102
|
35 // of segment a//b w.r.t. segment c
|
Chris@102
|
36 // This is calculated without (re)calculating the IP itself again and fully
|
Chris@102
|
37 // based on integer mathematics
|
Chris@102
|
38 template <typename T, typename Segment>
|
Chris@102
|
39 static inline T side_value(Segment const& a, Segment const& b,
|
Chris@102
|
40 Segment const& c)
|
Chris@102
|
41 {
|
Chris@102
|
42 // The first point of the three segments is reused several times
|
Chris@102
|
43 T const ax = get<0, 0>(a);
|
Chris@102
|
44 T const ay = get<0, 1>(a);
|
Chris@102
|
45 T const bx = get<0, 0>(b);
|
Chris@102
|
46 T const by = get<0, 1>(b);
|
Chris@102
|
47 T const cx = get<0, 0>(c);
|
Chris@102
|
48 T const cy = get<0, 1>(c);
|
Chris@102
|
49
|
Chris@102
|
50 T const dx_a = get<1, 0>(a) - ax;
|
Chris@102
|
51 T const dy_a = get<1, 1>(a) - ay;
|
Chris@102
|
52
|
Chris@102
|
53 T const dx_b = get<1, 0>(b) - bx;
|
Chris@102
|
54 T const dy_b = get<1, 1>(b) - by;
|
Chris@102
|
55
|
Chris@102
|
56 T const dx_c = get<1, 0>(c) - cx;
|
Chris@102
|
57 T const dy_c = get<1, 1>(c) - cy;
|
Chris@102
|
58
|
Chris@102
|
59 // Cramer's rule: d (see cart_intersect.hpp)
|
Chris@102
|
60 T const d = geometry::detail::determinant<T>
|
Chris@102
|
61 (
|
Chris@102
|
62 dx_a, dy_a,
|
Chris@102
|
63 dx_b, dy_b
|
Chris@102
|
64 );
|
Chris@102
|
65
|
Chris@102
|
66 T const zero = T();
|
Chris@102
|
67 if (d == zero)
|
Chris@102
|
68 {
|
Chris@102
|
69 // There is no IP of a//b, they are collinear or parallel
|
Chris@102
|
70 // We don't have to divide but we can already conclude the side-value
|
Chris@102
|
71 // is meaningless and the resulting determinant will be 0
|
Chris@102
|
72 return zero;
|
Chris@102
|
73 }
|
Chris@102
|
74
|
Chris@102
|
75 // Cramer's rule: da (see cart_intersect.hpp)
|
Chris@102
|
76 T const da = geometry::detail::determinant<T>
|
Chris@102
|
77 (
|
Chris@102
|
78 dx_b, dy_b,
|
Chris@102
|
79 ax - bx, ay - by
|
Chris@102
|
80 );
|
Chris@102
|
81
|
Chris@102
|
82 // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
|
Chris@102
|
83 // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
|
Chris@102
|
84 // We replace ipx by expression above and multiply each term by d
|
Chris@102
|
85 T const result = geometry::detail::determinant<T>
|
Chris@102
|
86 (
|
Chris@102
|
87 dx_c * d, dy_c * d,
|
Chris@102
|
88 d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
|
Chris@102
|
89 );
|
Chris@102
|
90
|
Chris@102
|
91 // Note: result / (d * d)
|
Chris@102
|
92 // is identical to the side_value of side_by_triangle
|
Chris@102
|
93 // Therefore, the sign is always the same as that result, and the
|
Chris@102
|
94 // resulting side (left,right,collinear) is the same
|
Chris@102
|
95
|
Chris@102
|
96 return result;
|
Chris@102
|
97
|
Chris@102
|
98 }
|
Chris@102
|
99
|
Chris@102
|
100 template <typename Segment>
|
Chris@102
|
101 static inline int apply(Segment const& a, Segment const& b, Segment const& c)
|
Chris@102
|
102 {
|
Chris@102
|
103 typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
|
Chris@102
|
104 coordinate_type const s = side_value<coordinate_type>(a, b, c);
|
Chris@102
|
105 coordinate_type const zero = coordinate_type();
|
Chris@102
|
106 return math::equals(s, zero) ? 0
|
Chris@102
|
107 : s > zero ? 1
|
Chris@102
|
108 : -1;
|
Chris@102
|
109 }
|
Chris@102
|
110
|
Chris@102
|
111 };
|
Chris@102
|
112
|
Chris@102
|
113
|
Chris@102
|
114 }} // namespace strategy::side
|
Chris@102
|
115
|
Chris@102
|
116 }} // namespace boost::geometry
|
Chris@102
|
117
|
Chris@102
|
118
|
Chris@102
|
119 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
|