1 /**
2 * \brief Various geometrical calculations.
3 */
5 #ifdef HAVE_CONFIG_H
6 # include <config.h>
7 #endif
8 #include <2geom/geom.h>
9 #include <2geom/point.h>
10 #include <algorithm>
11 #include <2geom/rect.h>
13 namespace Geom {
15 /**
16 * Finds the intersection of the two (infinite) lines
17 * defined by the points p such that dot(n0, p) == d0 and dot(n1, p) == d1.
18 *
19 * If the two lines intersect, then \a result becomes their point of
20 * intersection; otherwise, \a result remains unchanged.
21 *
22 * This function finds the intersection of the two lines (infinite)
23 * defined by n0.X = d0 and x1.X = d1. The algorithm is as follows:
24 * To compute the intersection point use kramer's rule:
25 * \verbatim
26 * convert lines to form
27 * ax + by = c
28 * dx + ey = f
29 *
30 * (
31 * e.g. a = (x2 - x1), b = (y2 - y1), c = (x2 - x1)*x1 + (y2 - y1)*y1
32 * )
33 *
34 * In our case we use:
35 * a = n0.x d = n1.x
36 * b = n0.y e = n1.y
37 * c = d0 f = d1
38 *
39 * so:
40 *
41 * adx + bdy = cd
42 * adx + aey = af
43 *
44 * bdy - aey = cd - af
45 * (bd - ae)y = cd - af
46 *
47 * y = (cd - af)/(bd - ae)
48 *
49 * repeat for x and you get:
50 *
51 * x = (fb - ce)/(bd - ae) \endverbatim
52 *
53 * If the denominator (bd-ae) is 0 then the lines are parallel, if the
54 * numerators are 0 then the lines coincide.
55 *
56 * \todo Why not use existing but outcommented code below
57 * (HAVE_NEW_INTERSECTOR_CODE)?
58 */
59 IntersectorKind
60 line_intersection(Geom::Point const &n0, double const d0,
61 Geom::Point const &n1, double const d1,
62 Geom::Point &result)
63 {
64 double denominator = dot(Geom::rot90(n0), n1);
65 double X = n1[Geom::Y] * d0 -
66 n0[Geom::Y] * d1;
67 /* X = (-d1, d0) dot (n0[Y], n1[Y]) */
69 if (denominator == 0) {
70 if ( X == 0 ) {
71 return coincident;
72 } else {
73 return parallel;
74 }
75 }
77 double Y = n0[Geom::X] * d1 -
78 n1[Geom::X] * d0;
80 result = Geom::Point(X, Y) / denominator;
82 return intersects;
83 }
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 the line segment from p00 to p01 intersects the
114 infinite line passing through p10 and p11. This doesn't find the
115 point of intersection, use the line_intersect function above,
116 or the segment_intersection interface below.
118 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
119 */
120 bool
121 line_segment_intersectp(Geom::Point const &p00, Geom::Point const &p01,
122 Geom::Point const &p10, Geom::Point const &p11)
123 {
124 if(p00 == p01) return false;
125 if(p10 == p11) return false;
127 return ((intersector_ccw(p00, p01, p10) * intersector_ccw(p00, p01, p11)) <= 0 );
128 }
131 /** Determine whether two line segments intersect. This doesn't find
132 the point of intersection, use the line_intersect function above,
133 or the segment_intersection interface below.
135 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
136 */
137 bool
138 segment_intersectp(Geom::Point const &p00, Geom::Point const &p01,
139 Geom::Point const &p10, Geom::Point const &p11)
140 {
141 if(p00 == p01) return false;
142 if(p10 == p11) return false;
144 /* true iff ( (the p1 segment straddles the p0 infinite line)
145 * and (the p0 segment straddles the p1 infinite line) ). */
146 return (line_segment_intersectp(p00, p01, p10, p11) &&
147 line_segment_intersectp(p10, p11, p00, p01));
148 }
150 /** Determine whether \& where a line segments intersects an (infinite) line.
152 If there is no intersection, then \a result remains unchanged.
154 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
155 **/
156 IntersectorKind
157 line_segment_intersect(Geom::Point const &p00, Geom::Point const &p01,
158 Geom::Point const &p10, Geom::Point const &p11,
159 Geom::Point &result)
160 {
161 if(line_segment_intersectp(p00, p01, p10, p11)) {
162 Geom::Point n0 = (p01 - p00).ccw();
163 double d0 = dot(n0,p00);
165 Geom::Point n1 = (p11 - p10).ccw();
166 double d1 = dot(n1,p10);
167 return line_intersection(n0, d0, n1, d1, result);
168 } else {
169 return no_intersection;
170 }
171 }
174 /** Determine whether \& where two line segments intersect.
176 If the two segments don't intersect, then \a result remains unchanged.
178 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
179 **/
180 IntersectorKind
181 segment_intersect(Geom::Point const &p00, Geom::Point const &p01,
182 Geom::Point const &p10, Geom::Point const &p11,
183 Geom::Point &result)
184 {
185 if(segment_intersectp(p00, p01, p10, p11)) {
186 Geom::Point n0 = (p01 - p00).ccw();
187 double d0 = dot(n0,p00);
189 Geom::Point n1 = (p11 - p10).ccw();
190 double d1 = dot(n1,p10);
191 return line_intersection(n0, d0, n1, d1, result);
192 } else {
193 return no_intersection;
194 }
195 }
197 /** Determine whether \& where two line segments intersect.
199 If the two segments don't intersect, then \a result remains unchanged.
201 \pre neither segment is zero-length; i.e. p00 != p01 and p10 != p11.
202 **/
203 IntersectorKind
204 line_twopoint_intersect(Geom::Point const &p00, Geom::Point const &p01,
205 Geom::Point const &p10, Geom::Point const &p11,
206 Geom::Point &result)
207 {
208 Geom::Point n0 = (p01 - p00).ccw();
209 double d0 = dot(n0,p00);
211 Geom::Point n1 = (p11 - p10).ccw();
212 double d1 = dot(n1,p10);
213 return line_intersection(n0, d0, n1, d1, result);
214 }
216 // this is used to compare points for std::sort below
217 static bool
218 is_less(Point const &A, Point const &B)
219 {
220 if (A[X] < B[X]) {
221 return true;
222 } else if (A[X] == B[X] && A[Y] < B[Y]) {
223 return true;
224 } else {
225 return false;
226 }
227 }
229 // TODO: this can doubtlessly be improved
230 static void
231 eliminate_duplicates_p(std::vector<Point> &pts)
232 {
233 unsigned int size = pts.size();
235 if (size < 2)
236 return;
238 if (size == 2) {
239 if (pts[0] == pts[1]) {
240 pts.pop_back();
241 }
242 } else {
243 std::sort(pts.begin(), pts.end(), &is_less);
244 if (size == 3) {
245 if (pts[0] == pts[1]) {
246 pts.erase(pts.begin());
247 } else if (pts[1] == pts[2]) {
248 pts.pop_back();
249 }
250 } else {
251 // we have size == 4
252 if (pts[2] == pts[3]) {
253 pts.pop_back();
254 }
255 if (pts[0] == pts[1]) {
256 pts.erase(pts.begin());
257 }
258 }
259 }
260 }
262 /** Determine whether \& where an (infinite) line intersects a rectangle.
263 *
264 * \a c0, \a c1 are diagonal corners of the rectangle and
265 * \a p1, \a p1 are distinct points on the line
266 *
267 * \return A list (possibly empty) of points of intersection. If two such points (say \a r0 and \a
268 * r1) then it is guaranteed that the order of \a r0, \a r1 along the line is the same as the that
269 * of \a c0, \a c1 (i.e., the vectors \a r1 - \a r0 and \a p1 - \a p0 point into the same
270 * direction).
271 */
272 std::vector<Geom::Point>
273 rect_line_intersect(Geom::Point const &c0, Geom::Point const &c1,
274 Geom::Point const &p0, Geom::Point const &p1)
275 {
276 using namespace Geom;
278 std::vector<Point> results;
280 Point A(c0);
281 Point C(c1);
283 Point B(A[X], C[Y]);
284 Point D(C[X], A[Y]);
286 Point res;
288 if (line_segment_intersect(p0, p1, A, B, res) == intersects) {
289 results.push_back(res);
290 }
291 if (line_segment_intersect(p0, p1, B, C, res) == intersects) {
292 results.push_back(res);
293 }
294 if (line_segment_intersect(p0, p1, C, D, res) == intersects) {
295 results.push_back(res);
296 }
297 if (line_segment_intersect(p0, p1, D, A, res) == intersects) {
298 results.push_back(res);
299 }
301 eliminate_duplicates_p(results);
303 if (results.size() == 2) {
304 // sort the results so that the order is the same as that of p0 and p1
305 Point dir1 (results[1] - results[0]);
306 Point dir2 (p1 - p0);
307 if (dot(dir1, dir2) < 0) {
308 std::swap(results[0], results[1]);
309 }
310 }
312 return results;
313 }
315 /** Determine whether \& where an (infinite) line intersects a rectangle.
316 *
317 * \a c0, \a c1 are diagonal corners of the rectangle and
318 * \a p1, \a p1 are distinct points on the line
319 *
320 * \return A list (possibly empty) of points of intersection. If two such points (say \a r0 and \a
321 * r1) then it is guaranteed that the order of \a r0, \a r1 along the line is the same as the that
322 * of \a c0, \a c1 (i.e., the vectors \a r1 - \a r0 and \a p1 - \a p0 point into the same
323 * direction).
324 */
325 boost::optional<LineSegment>
326 rect_line_intersect(Geom::Rect &r,
327 Geom::LineSegment ls)
328 {
329 std::vector<Point> results;
331 results = rect_line_intersect(r.min(), r.max(), ls[0], ls[1]);
332 if(results.size() == 2) {
333 return LineSegment(results[0], results[1]);
334 }
335 return boost::optional<LineSegment>();
336 }
338 boost::optional<LineSegment>
339 rect_line_intersect(Geom::Rect &r,
340 Geom::Line l)
341 {
342 return rect_line_intersect(r, l.segment(0, 1));
343 }
345 /**
346 * polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area of a polygon, given its
347 * vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It is assumed that the contour is closed, i.e., that
348 * the vertex following (x[n-1], y[n-1]) is (x[0], y[0]). The algebraic sign of the area is
349 * positive for counterclockwise ordering of vertices in x-y plane; otherwise negative.
351 * Returned values:
352 0 for normal execution;
353 1 if the polygon is degenerate (number of vertices < 3);
354 2 if area = 0 (and the centroid is undefined).
356 * for now we require the path to be a polyline and assume it is closed.
357 **/
359 int centroid(std::vector<Geom::Point> const &p, Geom::Point& centroid, double &area) {
360 const unsigned n = p.size();
361 if (n < 3)
362 return 1;
363 Geom::Point centroid_tmp(0,0);
364 double atmp = 0;
365 for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
366 const double ai = -cross(p[j], p[i]);
367 atmp += ai;
368 centroid_tmp += (p[j] + p[i])*ai; // first moment.
369 }
370 area = atmp / 2;
371 if (atmp != 0) {
372 centroid = centroid_tmp / (3 * atmp);
373 return 0;
374 }
375 return 2;
376 }
378 }
380 /*
381 Local Variables:
382 mode:c++
383 c-file-style:"stroustrup"
384 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
385 indent-tabs-mode:nil
386 fill-column:99
387 End:
388 */
389 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :