8845cefa89fe443ae321fe022a7ee4be8c22c574
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 "point.h"
12 namespace Geom {
14 /**
15 * Finds the intersection of the two (infinite) lines
16 * defined by the points p such that dot(n0, p) == d0 and dot(n1, p) == d1.
17 *
18 * If the two lines intersect, then \a result becomes their point of
19 * intersection; otherwise, \a result remains unchanged.
20 *
21 * This function finds the intersection of the two lines (infinite)
22 * defined by n0.X = d0 and x1.X = d1. The algorithm is as follows:
23 * To compute the intersection point use kramer's rule:
24 * \verbatim
25 * convert lines to form
26 * ax + by = c
27 * dx + ey = f
28 *
29 * (
30 * e.g. a = (x2 - x1), b = (y2 - y1), c = (x2 - x1)*x1 + (y2 - y1)*y1
31 * )
32 *
33 * In our case we use:
34 * a = n0.x d = n1.x
35 * b = n0.y e = n1.y
36 * c = d0 f = d1
37 *
38 * so:
39 *
40 * adx + bdy = cd
41 * adx + aey = af
42 *
43 * bdy - aey = cd - af
44 * (bd - ae)y = cd - af
45 *
46 * y = (cd - af)/(bd - ae)
47 *
48 * repeat for x and you get:
49 *
50 * x = (fb - ce)/(bd - ae) \endverbatim
51 *
52 * If the denominator (bd-ae) is 0 then the lines are parallel, if the
53 * numerators are then 0 then the lines coincide.
54 *
55 * \todo Why not use existing but outcommented code below
56 * (HAVE_NEW_INTERSECTOR_CODE)?
57 */
58 IntersectorKind
59 line_intersection(Geom::Point const &n0, double const d0,
60 Geom::Point const &n1, double const d1,
61 Geom::Point &result)
62 {
63 double denominator = dot(Geom::rot90(n0), n1);
64 double X = n1[Geom::Y] * d0 -
65 n0[Geom::Y] * d1;
66 /* X = (-d1, d0) dot (n0[Y], n1[Y]) */
68 if (denominator == 0) {
69 if ( X == 0 ) {
70 return coincident;
71 } else {
72 return parallel;
73 }
74 }
76 double Y = n0[Geom::X] * d1 -
77 n1[Geom::X] * d0;
79 result = Geom::Point(X, Y) / denominator;
81 return intersects;
82 }
87 /* ccw exists as a building block */
88 int
89 intersector_ccw(const Geom::Point& p0, const Geom::Point& p1,
90 const Geom::Point& p2)
91 /* Determine which way a set of three points winds. */
92 {
93 Geom::Point d1 = p1 - p0;
94 Geom::Point d2 = p2 - p0;
95 /* compare slopes but avoid division operation */
96 double c = dot(Geom::rot90(d1), d2);
97 if(c > 0)
98 return +1; // ccw - do these match def'n in header?
99 if(c < 0)
100 return -1; // cw
102 /* Colinear [or NaN]. Decide the order. */
103 if ( ( d1[0] * d2[0] < 0 ) ||
104 ( d1[1] * d2[1] < 0 ) ) {
105 return -1; // p2 < p0 < p1
106 } else if ( dot(d1,d1) < dot(d2,d2) ) {
107 return +1; // p0 <= p1 < p2
108 } else {
109 return 0; // p0 <= p2 <= p1
110 }
111 }
113 /** Determine whether two line segments intersect. This doesn't find
114 the point of intersection, use the line_intersect function above,
115 or the segment_intersection interface below.
117 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
118 */
119 static bool
120 segment_intersectp(Geom::Point const &p00, Geom::Point const &p01,
121 Geom::Point const &p10, Geom::Point const &p11)
122 {
123 if(p00 == p01) return false;
124 if(p10 == p11) return false;
126 /* true iff ( (the p1 segment straddles the p0 infinite line)
127 * and (the p0 segment straddles the p1 infinite line) ). */
128 return ((intersector_ccw(p00,p01, p10)
129 *intersector_ccw(p00, p01, p11)) <=0 )
130 &&
131 ((intersector_ccw(p10,p11, p00)
132 *intersector_ccw(p10, p11, p01)) <=0 );
133 }
136 /** Determine whether \& where two line segments intersect.
138 If the two segments don't intersect, then \a result remains unchanged.
140 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
141 **/
142 IntersectorKind
143 segment_intersect(Geom::Point const &p00, Geom::Point const &p01,
144 Geom::Point const &p10, Geom::Point const &p11,
145 Geom::Point &result)
146 {
147 if(segment_intersectp(p00, p01, p10, p11)) {
148 Geom::Point n0 = (p01 - p00).ccw();
149 double d0 = dot(n0,p00);
151 Geom::Point n1 = (p11 - p10).ccw();
152 double d1 = dot(n1,p10);
153 return line_intersection(n0, d0, n1, d1, result);
154 } else {
155 return no_intersection;
156 }
157 }
159 /** Determine whether \& where two line segments intersect.
161 If the two segments don't intersect, then \a result remains unchanged.
163 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
164 **/
165 IntersectorKind
166 line_twopoint_intersect(Geom::Point const &p00, Geom::Point const &p01,
167 Geom::Point const &p10, Geom::Point const &p11,
168 Geom::Point &result)
169 {
170 Geom::Point n0 = (p01 - p00).ccw();
171 double d0 = dot(n0,p00);
173 Geom::Point n1 = (p11 - p10).ccw();
174 double d1 = dot(n1,p10);
175 return line_intersection(n0, d0, n1, d1, result);
176 }
178 /**
179 * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its
180 * vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It is assumed that the contour is closed, i.e., that
181 * the vertex following (x[n-1], y[n-1]) is (x[0], y[0]). The algebraic sign of the area is
182 * positive for counterclockwise ordering of vertices in x-y plane; otherwise negative.
184 * Returned values:
185 0 for normal execution;
186 1 if the polygon is degenerate (number of vertices < 3);
187 2 if area = 0 (and the centroid is undefined).
189 * for now we require the path to be a polyline and assume it is closed.
190 **/
192 int centroid(std::vector<Geom::Point> p, Geom::Point& centroid, double &area) {
193 const unsigned n = p.size();
194 if (n < 3)
195 return 1;
196 Geom::Point centroid_tmp(0,0);
197 double atmp = 0;
198 for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
199 const double ai = -cross(p[j], p[i]);
200 atmp += ai;
201 centroid_tmp += (p[j] + p[i])*ai; // first moment.
202 }
203 area = atmp / 2;
204 if (atmp != 0) {
205 centroid = centroid_tmp / (3 * atmp);
206 return 0;
207 }
208 return 2;
209 }
211 }
213 /*
214 Local Variables:
215 mode:c++
216 c-file-style:"stroustrup"
217 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
218 indent-tabs-mode:nil
219 fill-column:99
220 End:
221 */
222 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :