1 /**
2 * \file src/geom.cpp
3 * \brief Various geometrical calculations.
4 */
6 #ifdef HAVE_CONFIG_H
7 # include <config.h>
8 #endif
9 #include "geom.h"
10 #include <libnr/nr-point-fns.h>
12 /**
13 * Finds the intersection of the two (infinite) lines
14 * defined by the points p such that dot(n0, p) == d0 and dot(n1, p) == d1.
15 *
16 * If the two lines intersect, then \a result becomes their point of
17 * intersection; otherwise, \a result remains unchanged.
18 *
19 * This function finds the intersection of the two lines (infinite)
20 * defined by n0.X = d0 and x1.X = d1. The algorithm is as follows:
21 * To compute the intersection point use kramer's rule:
22 * \verbatim
23 * convert lines to form
24 * ax + by = c
25 * dx + ey = f
26 *
27 * (
28 * e.g. a = (x2 - x1), b = (y2 - y1), c = (x2 - x1)*x1 + (y2 - y1)*y1
29 * )
30 *
31 * In our case we use:
32 * a = n0.x d = n1.x
33 * b = n0.y e = n1.y
34 * c = d0 f = d1
35 *
36 * so:
37 *
38 * adx + bdy = cd
39 * adx + aey = af
40 *
41 * bdy - aey = cd - af
42 * (bd - ae)y = cd - af
43 *
44 * y = (cd - af)/(bd - ae)
45 *
46 * repeat for x and you get:
47 *
48 * x = (fb - ce)/(bd - ae) \endverbatim
49 *
50 * If the denominator (bd-ae) is 0 then the lines are parallel, if the
51 * numerators are then 0 then the lines coincide.
52 *
53 * \todo Why not use existing but outcommented code below
54 * (HAVE_NEW_INTERSECTOR_CODE)?
55 */
56 IntersectorKind intersector_line_intersection(NR::Point const &n0, double const d0,
57 NR::Point const &n1, double const d1,
58 NR::Point &result)
59 {
60 double denominator = dot(NR::rot90(n0), n1);
61 double X = n1[NR::Y] * d0 -
62 n0[NR::Y] * d1;
63 /* X = (-d1, d0) dot (n0[Y], n1[Y]) */
65 if (denominator == 0) {
66 if ( X == 0 ) {
67 return COINCIDENT;
68 } else {
69 return PARALLEL;
70 }
71 }
73 double Y = n0[NR::X] * d1 -
74 n1[NR::X] * d0;
76 result = NR::Point(X, Y) / denominator;
78 return INTERSECTS;
79 }
84 /*
85 New code which we are not yet using
86 */
87 #ifdef HAVE_NEW_INTERSECTOR_CODE
90 /* ccw exists as a building block */
91 static int
92 sp_intersector_ccw(const NR::Point p0, const NR::Point p1, const NR::Point p2)
93 /* Determine which way a set of three points winds. */
94 {
95 NR::Point d1 = p1 - p0;
96 NR::Point d2 = p2 - p0;
97 /* compare slopes but avoid division operation */
98 double c = dot(NR::rot90(d1), d2);
99 if(c > 0)
100 return +1; // ccw - do these match def'n in header?
101 if(c < 0)
102 return -1; // cw
104 /* Colinear [or NaN]. Decide the order. */
105 if ( ( d1[0] * d2[0] < 0 ) ||
106 ( d1[1] * d2[1] < 0 ) ) {
107 return -1; // p2 < p0 < p1
108 } else if ( dot(d1,d1) < dot(d2,d2) ) {
109 return +1; // p0 <= p1 < p2
110 } else {
111 return 0; // p0 <= p2 <= p1
112 }
113 }
115 /** Determine whether two line segments intersect. This doesn't find
116 the point of intersection, use the line_intersect function above,
117 or the segment_intersection interface below.
119 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
120 */
121 static bool
122 sp_intersector_segment_intersectp(NR::Point const &p00, NR::Point const &p01,
123 NR::Point const &p10, NR::Point const &p11)
124 {
125 g_return_val_if_fail(p00 != p01, false);
126 g_return_val_if_fail(p10 != p11, false);
128 /* true iff ( (the p1 segment straddles the p0 infinite line)
129 * and (the p0 segment straddles the p1 infinite line) ). */
130 return ((sp_intersector_ccw(p00,p01, p10)
131 *sp_intersector_ccw(p00, p01, p11)) <=0 )
132 &&
133 ((sp_intersector_ccw(p10,p11, p00)
134 *sp_intersector_ccw(p10, p11, p01)) <=0 );
135 }
138 /** Determine whether \& where two line segments intersect.
140 If the two segments don't intersect, then \a result remains unchanged.
142 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
143 **/
144 static sp_intersector_kind
145 sp_intersector_segment_intersect(NR::Point const &p00, NR::Point const &p01,
146 NR::Point const &p10, NR::Point const &p11,
147 NR::Point &result)
148 {
149 if(sp_intersector_segment_intersectp(p00, p01, p10, p11)) {
150 NR::Point n0 = (p00 - p01).ccw();
151 double d0 = dot(n0,p00);
153 NR::Point n1 = (p10 - p11).ccw();
154 double d1 = dot(n1,p10);
155 return sp_intersector_line_intersection(n0, d0, n1, d1, result);
156 } else {
157 return no_intersection;
158 }
159 }
161 #endif /* end yet-unused HAVE_NEW_INTERSECTOR_CODE code */
163 /*
164 Local Variables:
165 mode:c++
166 c-file-style:"stroustrup"
167 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
168 indent-tabs-mode:nil
169 fill-column:99
170 End:
171 */
172 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :