1
// Generic Geometry Library
3
// Copyright Barend Gehrels 1995-2009, Geodan Holding B.V. Amsterdam, the Netherlands.
4
// Copyright Bruno Lalande 2008, 2009
5
// Use, modification and distribution is subject to the Boost Software License,
6
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
7
// http://www.boost.org/LICENSE_1_0.txt)
9
#ifndef GGL_STRATEGY_CARTESIAN_INTERSECTION_HPP
10
#define GGL_STRATEGY_CARTESIAN_INTERSECTION_HPP
14
#include <ggl/core/concepts/point_concept.hpp>
15
#include <ggl/core/concepts/segment_concept.hpp>
21
namespace strategy { namespace intersection {
24
#ifndef DOXYGEN_NO_DETAIL
28
template <typename S, size_t D>
29
struct segment_interval
32
static inline void arrange(S const& s, T& s_1, T& s_2, bool& swapped)
49
\see http://mathworld.wolfram.com/Line-LineIntersection.html
52
struct relate_cartesian_segments
54
typedef typename F::return_type RETURN_TYPE;
55
typedef typename F::segment_type1 S1;
56
typedef typename F::segment_type2 S2;
58
typedef typename point_type<S1>::type P;
59
BOOST_CONCEPT_ASSERT( (concept::Point<P>) );
61
BOOST_CONCEPT_ASSERT( (concept::ConstSegment<S1>) );
62
BOOST_CONCEPT_ASSERT( (concept::ConstSegment<S2>) );
63
typedef typename select_coordinate_type<S1, S2>::type T;
65
/// Relate segments a and b
66
static inline RETURN_TYPE relate(S1 const& a, S2 const& b)
68
T dx_a = get<1, 0>(a) - get<0, 0>(a); // distance in x-dir
69
T dx_b = get<1, 0>(b) - get<0, 0>(b);
70
T dy_a = get<1, 1>(a) - get<0, 1>(a); // distance in y-dir
71
T dy_b = get<1, 1>(b) - get<0, 1>(b);
72
return relate(a, b, dx_a, dy_a, dx_b, dy_b);
76
/// Relate segments a and b using precalculated differences. This can save two or four substractions in many cases
77
static inline RETURN_TYPE relate(S1 const& a, S2 const& b,
78
T const& dx_a, T const& dy_a, T const& dx_b, T const& dy_b)
80
T wx = get<0, 0>(a) - get<0, 0>(b);
81
T wy = get<0, 1>(a) - get<0, 1>(b);
83
// Calculate determinants - Cramers rule
84
T d = (dy_b * dx_a) - (dx_b * dy_a);
85
T da = (dx_b * wy) - (dy_b * wx);
86
T db = (dx_a * wy) - (dy_a * wx);
88
if(! math::equals(d, 0))
90
// Determinant d is nonzero. Rays do intersect. This is the normal case.
91
// ra/rb: ratio 0-1 where intersection divides A/B
94
// The ra/rb calculation can probably also be avoided, only necessary to calculate the points themselves
95
// On segment: 0 <= r <= 1
96
// where: r==0 and r==1 are special cases
97
// --> r=0 if and only if da==0 (d != 0) [this is also checked below]
99
// --> da/d > 0 if da==positive and d==positive OR da==neg and d==neg
100
// --> da/d < 1 if (both positive) da < d or (negative) da > d, e.g. -2 > -4
101
// --> it would save a division but makes it complexer
103
double ra = double(da) / double(d);
104
double rb = double(db) / double(d);
106
return F::rays_intersect(ra >= 0 && ra <= 1 && rb >= 0 && rb <= 1, ra, rb,
107
dx_a, dy_a, dx_b, dy_b, wx, wy, a, b);
110
if(math::equals(da, 0) && math::equals(db, 0))
112
// Both da & db are zero. Segments are collinear. We'll find out how.
113
return relate_collinear(a, b, dx_a, dy_a, dx_b, dy_b);
116
// Segments are parallel (might still be opposite but this is probably never interesting)
117
return F::parallel();
123
/// Relate segments known collinear
124
static inline RETURN_TYPE relate_collinear(S1 const& a, S2 const& b,
125
T const& dx_a, T const& dy_a,
126
T const& dx_b, T const& dy_b)
128
// All ca. 200 lines are about collinear rays
129
// The intersections, if any, are always boundary points of the segments. No need to calculate anything.
130
// However we want to find out HOW they intersect, there are many cases.
131
// Most sources only provide the intersection (above) or that there is a collinearity (but not the points)
132
// or some spare sources give the intersection points (calculated) but not how they align.
133
// This source tries to give everything and still be efficient.
134
// It is therefore (and because of the extensive clarification comments) rather long...
136
// \see http://mpa.itc.it/radim/g50history/CMP/4.2.1-CERL-beta-libes/file475.txt
137
// \see http://docs.codehaus.org/display/GEOTDOC/Point+Set+Theory+and+the+DE-9IM+Matrix
138
// \see http://mathworld.wolfram.com/Line-LineIntersection.html
140
// Because of collinearity the case is now one-dimensional and can be checked using intervals
141
// We arrange either horizontally or vertically
142
// We get then two intervals:
143
// a_1-------------a_2 where a_1 < a_2
144
// b_1-------------b_2 where b_1 < b_2
145
// In all figures below a_1/a_2 denotes arranged intervals, a1-a2 or a2-a1 are still unarranged
146
T a_1, a_2, b_1, b_2;
147
bool a_swapped = false, b_swapped = false;
148
if (math::equals(dx_b, 0))
150
// Vertical -> Check y-direction
151
detail::segment_interval<S1, 1>::arrange(a, a_1, a_2, a_swapped);
152
detail::segment_interval<S1, 1>::arrange(b, b_1, b_2, b_swapped);
157
detail::segment_interval<S1, 0>::arrange(a, a_1, a_2, a_swapped);
158
detail::segment_interval<S1, 0>::arrange(b, b_1, b_2, b_swapped);
161
// First handle "disjoint", probably common case.
162
// 2 cases: a_1----------a_2 b_1-------b_2 or B left of A
163
if (a_2 < b_1 || a_1 > b_2)
165
return F::collinear_disjoint();
168
// Then handle "equal", in polygon neighbourhood comparisons also a common case
170
// Check if segments are equal...
171
bool a1_eq_b1 = math::equals(get<0, 0>(a), get<0, 0>(b))
172
&& math::equals(get<0, 1>(a), get<0, 1>(b));
173
bool a2_eq_b2 = math::equals(get<1, 0>(a), get<1, 0>(b))
174
&& math::equals(get<1, 1>(a), get<1, 1>(b));
175
if (a1_eq_b1 && a2_eq_b2)
177
return F::segment_equal(a, false);
180
// ... or opposite equal
181
bool a1_eq_b2 = math::equals(get<0, 0>(a), get<1, 0>(b))
182
&& math::equals(get<0, 1>(a), get<1, 1>(b));
183
bool a2_eq_b1 = math::equals(get<1, 0>(a), get<0, 0>(b))
184
&& math::equals(get<1, 1>(a), get<0, 1>(b));
185
if (a1_eq_b2 && a2_eq_b1)
187
return F::segment_equal(a, true);
191
// Degenerate cases: segments of single point, lying on other segment, non disjoint
192
if (math::equals(dx_a, 0) && math::equals(dy_a, 0))
194
return F::degenerate(a, true);
196
if (math::equals(dx_b, 0) && math::equals(dy_b, 0))
198
return F::degenerate(b, false);
202
// The rest below will return one or two intersections.
203
// The delegated class can decide which is the intersection point, or two, build the Intersection Matrix (IM)
204
// For IM it is important to know which relates to which. So this information is given,
205
// without performance penalties to intersection calculation
207
bool has_common_points = a1_eq_b1 || a1_eq_b2 || a2_eq_b1 || a2_eq_b2;
210
// "Touch" -> one intersection point -> one but not two common points
211
// --------> A (or B)
212
// <---------- B (or A)
213
// a_2==b_1 (b_2==a_1 or a_2==b1)
215
// The check a_2/b_1 is necessary because it excludes cases like
218
// ... which are handled lateron
220
// Corresponds to 4 cases, of which the equal points are determined above
221
// 1: a1---->a2 b1--->b2
222
// 2: a2<----a1 b2<---b1
223
// 3: a1---->a2 b2<---b1
224
// 4: a2<----a1 b1--->b2
225
// Where the arranged forms have two forms:
226
// a_1-----a_2/b_1-------b_2 or reverse (B left of A)
227
if (has_common_points && (math::equals(a_2, b_1) || math::equals(b_2, a_1)))
229
if (a2_eq_b1) return F::collinear_touch(get<1, 0>(a), get<1, 1>(a), false);
230
if (a2_eq_b2) return F::collinear_touch(get<1, 0>(a), get<1, 1>(a), true);
231
if (a1_eq_b2) return F::collinear_touch(get<0, 0>(a), get<0, 1>(a), false);
232
if (a1_eq_b1) return F::collinear_touch(get<0, 0>(a), get<0, 1>(a), true);
236
// "Touch/within" -> there are common points and also an intersection of interiors:
237
// Corresponds to many cases:
238
// Case 1: a1------->a2 Case 5: a1-->a2 Case 9: a1--->a2
239
// b1--->b2 b1------->b2 b1---------b2
240
// Case 2: a2<-------a1
241
// b1----b2 Et cetera
242
// Case 3: a1------->a2
244
// Case 4: a2<-------a1
247
// For case 1-4: a_1 < (b_1 or b_2) < a_2, two intersections are equal to segment B
248
// For case 5-8: b_1 < (a_1 or a_2) < b_2, two intersections are equal to segment A
249
if (has_common_points)
251
bool a_in_b = (b_1 < a_1 && a_1 < b_2) || (b_1 < a_2 && a_2 < b_2);
252
if (a2_eq_b2) return F::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, false);
253
if (a1_eq_b2) return F::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, true);
254
if (a2_eq_b1) return F::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, true);
255
if (a1_eq_b1) return F::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, false);
258
bool opposite = a_swapped ^ b_swapped;
261
// "Inside", a completely within b or b completely within a
264
// a_1---a_2 -> take A's points as intersection points
265
// b_1------------b_2
267
// a_1------------a_2
268
// b_1---b_2 -> take B's points
269
if (a_1 > b_1 && a_2 < b_2)
272
return F::collinear_a_in_b(a, opposite);
274
if (b_1 > a_1 && b_2 < a_2)
277
return F::collinear_b_in_a(b, opposite);
281
// Now that all cases with equal,touch,inside,disjoint,degenerate are handled the only thing left is an overlap
282
// Case 1: a1--------->a2 a_1---------a_2
283
// b1----->b2 b_1------b_2
284
// Case 2: a2<---------a1 a_1---------a_2 a_swapped
285
// b1----->b2 b_1------b_2
286
// Case 3: a1--------->a2 a_1---------a_2
287
// b2<-----b1 b_1------b_2 b_swapped
288
// Case 4: a2<-------->a1 a_1---------a_2 a_swapped
289
// b2<-----b1 b_1------b_2 b_swapped
291
// Case 5: a1--------->a2 a_1---------a_2
292
// b1----->b2 b1--------b2
293
// Case 6: a2<---------a1
295
// Case 7: a1--------->a2
297
// Case 8: a2<-------->a1
300
if (a_1 < b_1 && b_1 < a_2)
304
! a_swapped && ! b_swapped ? F::collinear_overlaps(get<1, 0>(a), get<1, 1>(a), get<0, 0>(b), get<0, 1>(b), opposite)
305
: a_swapped ? F::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<0, 0>(b), get<0, 1>(b), opposite)
306
: b_swapped ? F::collinear_overlaps(get<1, 0>(a), get<1, 1>(a), get<1, 0>(b), get<1, 1>(b), opposite)
307
: F::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<1, 0>(b), get<1, 1>(b), opposite)
310
if (b_1 < a_1 && a_1 < b_2)
314
! a_swapped && ! b_swapped ? F::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<1, 0>(b), get<1, 1>(b), opposite)
315
: a_swapped ? F::collinear_overlaps(get<1, 0>(a), get<1, 1>(a), get<1, 0>(b), get<1, 1>(b), opposite)
316
: b_swapped ? F::collinear_overlaps(get<0, 0>(a), get<0, 1>(a), get<0, 0>(b), get<0, 1>(b), opposite)
317
: F::collinear_overlaps(get<1, 0>(a), get<1, 1>(a), get<0, 0>(b), get<0, 1>(b), opposite)
321
// Nothing should goes through. If any we have made an error
322
// TODO: proper exception
323
throw std::string("Internal error, unexpected case in relate_segment");
328
}} // namespace strategy::intersection
335
#endif // GGL_STRATEGY_CARTESIAN_INTERSECTION_HPP
added
removed
include/ggl/strategies/cartesian/cart_intersect.hpp
