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 Cramer's rule:
22 * (see http://en.wikipedia.org/wiki/Cramer%27s_rule)
23 * \verbatim
24 * convert lines to form
25 * ax + by = c
26 * dx + ey = f
27 *
28 * (
29 * e.g. a = (x2 - x1), b = (y2 - y1), c = (x2 - x1)*x1 + (y2 - y1)*y1
30 * )
31 *
32 * In our case we use:
33 * a = n0.x d = n1.x
34 * b = n0.y e = n1.y
35 * c = d0 f = d1
36 *
37 * so:
38 *
39 * adx + bdy = cd
40 * adx + aey = af
41 *
42 * bdy - aey = cd - af
43 * (bd - ae)y = cd - af
44 *
45 * y = (cd - af)/(bd - ae)
46 *
47 * repeat for x and you get:
48 *
49 * x = (fb - ce)/(bd - ae) \endverbatim
50 *
51 * If the denominator (bd-ae) is 0 then the lines are parallel, if the
52 * numerators are then 0 then the lines coincide.
53 *
54 * \todo Why not use existing but outcommented code below
55 * (HAVE_NEW_INTERSECTOR_CODE)?
56 */
57 IntersectorKind intersector_line_intersection(NR::Point const &n0, double const d0,
58 NR::Point const &n1, double const d1,
59 NR::Point &result)
60 {
61 double denominator = dot(NR::rot90(n0), n1);
62 double X = n1[NR::Y] * d0 -
63 n0[NR::Y] * d1;
64 /* X = (-d1, d0) dot (n0[Y], n1[Y]) */
66 if (denominator == 0) {
67 if ( X == 0 ) {
68 return COINCIDENT;
69 } else {
70 return PARALLEL;
71 }
72 }
74 double Y = n0[NR::X] * d1 -
75 n1[NR::X] * d0;
77 result = NR::Point(X, Y) / denominator;
79 return INTERSECTS;
80 }
85 /*
86 New code which we are not yet using
87 */
88 #ifdef HAVE_NEW_INTERSECTOR_CODE
91 /* ccw exists as a building block */
92 static int
93 sp_intersector_ccw(const NR::Point p0, const NR::Point p1, const NR::Point p2)
94 /* Determine which way a set of three points winds. */
95 {
96 NR::Point d1 = p1 - p0;
97 NR::Point d2 = p2 - p0;
98 /* compare slopes but avoid division operation */
99 double c = dot(NR::rot90(d1), d2);
100 if(c > 0)
101 return +1; // ccw - do these match def'n in header?
102 if(c < 0)
103 return -1; // cw
105 /* Colinear [or NaN]. Decide the order. */
106 if ( ( d1[0] * d2[0] < 0 ) ||
107 ( d1[1] * d2[1] < 0 ) ) {
108 return -1; // p2 < p0 < p1
109 } else if ( dot(d1,d1) < dot(d2,d2) ) {
110 return +1; // p0 <= p1 < p2
111 } else {
112 return 0; // p0 <= p2 <= p1
113 }
114 }
116 /** Determine whether two line segments intersect. This doesn't find
117 the point of intersection, use the line_intersect function above,
118 or the segment_intersection interface below.
120 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
121 */
122 static bool
123 sp_intersector_segment_intersectp(NR::Point const &p00, NR::Point const &p01,
124 NR::Point const &p10, NR::Point const &p11)
125 {
126 g_return_val_if_fail(p00 != p01, false);
127 g_return_val_if_fail(p10 != p11, false);
129 /* true iff ( (the p1 segment straddles the p0 infinite line)
130 * and (the p0 segment straddles the p1 infinite line) ). */
131 return ((sp_intersector_ccw(p00,p01, p10)
132 *sp_intersector_ccw(p00, p01, p11)) <=0 )
133 &&
134 ((sp_intersector_ccw(p10,p11, p00)
135 *sp_intersector_ccw(p10, p11, p01)) <=0 );
136 }
139 /** Determine whether \& where two line segments intersect.
141 If the two segments don't intersect, then \a result remains unchanged.
143 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
144 **/
145 static sp_intersector_kind
146 sp_intersector_segment_intersect(NR::Point const &p00, NR::Point const &p01,
147 NR::Point const &p10, NR::Point const &p11,
148 NR::Point &result)
149 {
150 if(sp_intersector_segment_intersectp(p00, p01, p10, p11)) {
151 NR::Point n0 = (p00 - p01).ccw();
152 double d0 = dot(n0,p00);
154 NR::Point n1 = (p10 - p11).ccw();
155 double d1 = dot(n1,p10);
156 return sp_intersector_line_intersection(n0, d0, n1, d1, result);
157 } else {
158 return no_intersection;
159 }
160 }
162 #endif /* end yet-unused HAVE_NEW_INTERSECTOR_CODE code */
164 /*
165 Local Variables:
166 mode:c++
167 c-file-style:"stroustrup"
168 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
169 indent-tabs-mode:nil
170 fill-column:99
171 End:
172 */
173 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :