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