1 #include "basic-intersection.h"
2 #include "exception.h"
4 unsigned intersect_steps = 0;
6 using std::vector;
7 namespace Geom {
9 class OldBezier {
10 public:
11 std::vector<Geom::Point> p;
12 OldBezier() {
13 }
14 void split(double t, OldBezier &a, OldBezier &b) const;
16 ~OldBezier() {}
18 void bounds(double &minax, double &maxax,
19 double &minay, double &maxay) {
20 // Compute bounding box for a
21 minax = p[0][X]; // These are the most likely to be extremal
22 maxax = p.back()[X];
23 if( minax > maxax )
24 std::swap(minax, maxax);
25 for(unsigned i = 1; i < p.size()-1; i++) {
26 if( p[i][X] < minax )
27 minax = p[i][X];
28 else if( p[i][X] > maxax )
29 maxax = p[i][X];
30 }
32 minay = p[0][Y]; // These are the most likely to be extremal
33 maxay = p.back()[Y];
34 if( minay > maxay )
35 std::swap(minay, maxay);
36 for(unsigned i = 1; i < p.size()-1; i++) {
37 if( p[i][Y] < minay )
38 minay = p[i][Y];
39 else if( p[i][Y] > maxay )
40 maxay = p[i][Y];
41 }
43 }
45 };
47 static std::vector<std::pair<double, double> >
48 find_intersections( OldBezier a, OldBezier b);
50 static std::vector<std::pair<double, double> >
51 find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A);
53 std::vector<std::pair<double, double> >
54 find_intersections( vector<Geom::Point> const & A,
55 vector<Geom::Point> const & B) {
56 OldBezier a, b;
57 a.p = A;
58 b.p = B;
59 return find_intersections(a,b);
60 }
62 std::vector<std::pair<double, double> >
63 find_self_intersections(OldBezier const &Sb) {
64 throwNotImplemented();
65 }
67 std::vector<std::pair<double, double> >
68 find_self_intersections(D2<SBasis> const & A) {
69 OldBezier Sb;
70 Sb.p = sbasis_to_bezier(A);
71 return find_self_intersections(Sb, A);
72 }
75 static std::vector<std::pair<double, double> >
76 find_self_intersections(OldBezier const &Sb, D2<SBasis> const & A) {
79 vector<double> dr = roots(derivative(A[X]));
80 {
81 vector<double> dyr = roots(derivative(A[Y]));
82 dr.insert(dr.begin(), dyr.begin(), dyr.end());
83 }
84 dr.push_back(0);
85 dr.push_back(1);
86 // We want to be sure that we have no empty segments
87 sort(dr.begin(), dr.end());
88 unique(dr.begin(), dr.end());
90 std::vector<std::pair<double, double> > all_si;
92 vector<OldBezier> pieces;
93 {
94 OldBezier in = Sb, l, r;
95 for(unsigned i = 0; i < dr.size()-1; i++) {
96 in.split((dr[i+1]-dr[i]) / (1 - dr[i]), l, r);
97 pieces.push_back(l);
98 in = r;
99 }
100 }
101 for(unsigned i = 0; i < dr.size()-1; i++) {
102 for(unsigned j = i+1; j < dr.size()-1; j++) {
103 std::vector<std::pair<double, double> > section =
104 find_intersections( pieces[i], pieces[j]);
105 for(unsigned k = 0; k < section.size(); k++) {
106 double l = section[k].first;
107 double r = section[k].second;
108 // XXX: This condition will prune out false positives, but it might create some false negatives. Todo: Confirm it is correct.
109 if(j == i+1)
110 if((l == 1) && (r == 0))
111 continue;
112 all_si.push_back(std::make_pair((1-l)*dr[i] + l*dr[i+1],
113 (1-r)*dr[j] + r*dr[j+1]));
114 }
115 }
116 }
118 // Because i is in order, all_si should be roughly already in order?
119 //sort(all_si.begin(), all_si.end());
120 //unique(all_si.begin(), all_si.end());
122 return all_si;
123 }
125 /* The value of 1.0 / (1L<<14) is enough for most applications */
126 const double INV_EPS = (1L<<14);
128 /*
129 * split the curve at the midpoint, returning an array with the two parts
130 * Temporary storage is minimized by using part of the storage for the result
131 * to hold an intermediate value until it is no longer needed.
132 */
133 void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
134 const unsigned sz = p.size();
135 Geom::Point Vtemp[sz][sz];
137 /* Copy control points */
138 std::copy(p.begin(), p.end(), Vtemp[0]);
140 /* Triangle computation */
141 for (unsigned i = 1; i < sz; i++) {
142 for (unsigned j = 0; j < sz - i; j++) {
143 Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
144 }
145 }
147 left.p.resize(sz);
148 right.p.resize(sz);
149 for (unsigned j = 0; j < sz; j++)
150 left.p[j] = Vtemp[j][0];
151 for (unsigned j = 0; j < sz; j++)
152 right.p[j] = Vtemp[sz-1-j][j];
153 }
156 /*
157 * Test the bounding boxes of two OldBezier curves for interference.
158 * Several observations:
159 * First, it is cheaper to compute the bounding box of the second curve
160 * and test its bounding box for interference than to use a more direct
161 * approach of comparing all control points of the second curve with
162 * the various edges of the bounding box of the first curve to test
163 * for interference.
164 * Second, after a few subdivisions it is highly probable that two corners
165 * of the bounding box of a given Bezier curve are the first and last
166 * control point. Once this happens once, it happens for all subsequent
167 * subcurves. It might be worth putting in a test and then short-circuit
168 * code for further subdivision levels.
169 * Third, in the final comparison (the interference test) the comparisons
170 * should both permit equality. We want to find intersections even if they
171 * occur at the ends of segments.
172 * Finally, there are tighter bounding boxes that can be derived. It isn't
173 * clear whether the higher probability of rejection (and hence fewer
174 * subdivisions and tests) is worth the extra work.
175 */
177 bool intersect_BB( OldBezier a, OldBezier b ) {
178 double minax, maxax, minay, maxay;
179 a.bounds(minax, maxax, minay, maxay);
180 double minbx, maxbx, minby, maxby;
181 b.bounds(minbx, maxbx, minby, maxby);
182 // Test bounding box of b against bounding box of a
183 // Not >= : need boundary case
184 return not( ( minax > maxbx ) || ( minay > maxby )
185 || ( minbx > maxax ) || ( minby > maxay ) );
186 }
188 /*
189 * Recursively intersect two curves keeping track of their real parameters
190 * and depths of intersection.
191 * The results are returned in a 2-D array of doubles indicating the parameters
192 * for which intersections are found. The parameters are in the order the
193 * intersections were found, which is probably not in sorted order.
194 * When an intersection is found, the parameter value for each of the two
195 * is stored in the index elements array, and the index is incremented.
196 *
197 * If either of the curves has subdivisions left before it is straight
198 * (depth > 0)
199 * that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
200 * the depth(s) is (are) decremented, and the parameter value(s) corresponding
201 * to the midpoints(s) is (are) computed.
202 * Then each of the subcurves of one curve is intersected with each of the
203 * subcurves of the other curve, first by testing the bounding boxes for
204 * interference. If there is any bounding box interference, the corresponding
205 * subcurves are recursively intersected.
206 *
207 * If neither curve has subdivisions left, the line segments from the first
208 * to last control point of each segment are intersected. (Actually the
209 * only the parameter value corresponding to the intersection point is found).
210 *
211 * The apriori flatness test is probably more efficient than testing at each
212 * level of recursion, although a test after three or four levels would
213 * probably be worthwhile, since many curves become flat faster than their
214 * asymptotic rate for the first few levels of recursion.
215 *
216 * The bounding box test fails much more frequently than it succeeds, providing
217 * substantial pruning of the search space.
218 *
219 * Each (sub)curve is subdivided only once, hence it is not possible that for
220 * one final line intersection test the subdivision was at one level, while
221 * for another final line intersection test the subdivision (of the same curve)
222 * was at another. Since the line segments share endpoints, the intersection
223 * is robust: a near-tangential intersection will yield zero or two
224 * intersections.
225 */
226 void recursively_intersect( OldBezier a, double t0, double t1, int deptha,
227 OldBezier b, double u0, double u1, int depthb,
228 std::vector<std::pair<double, double> > ¶meters)
229 {
230 intersect_steps ++;
231 if( deptha > 0 )
232 {
233 OldBezier A[2];
234 a.split(0.5, A[0], A[1]);
235 double tmid = (t0+t1)*0.5;
236 deptha--;
237 if( depthb > 0 )
238 {
239 OldBezier B[2];
240 b.split(0.5, B[0], B[1]);
241 double umid = (u0+u1)*0.5;
242 depthb--;
243 if( intersect_BB( A[0], B[0] ) )
244 recursively_intersect( A[0], t0, tmid, deptha,
245 B[0], u0, umid, depthb,
246 parameters );
247 if( intersect_BB( A[1], B[0] ) )
248 recursively_intersect( A[1], tmid, t1, deptha,
249 B[0], u0, umid, depthb,
250 parameters );
251 if( intersect_BB( A[0], B[1] ) )
252 recursively_intersect( A[0], t0, tmid, deptha,
253 B[1], umid, u1, depthb,
254 parameters );
255 if( intersect_BB( A[1], B[1] ) )
256 recursively_intersect( A[1], tmid, t1, deptha,
257 B[1], umid, u1, depthb,
258 parameters );
259 }
260 else
261 {
262 if( intersect_BB( A[0], b ) )
263 recursively_intersect( A[0], t0, tmid, deptha,
264 b, u0, u1, depthb,
265 parameters );
266 if( intersect_BB( A[1], b ) )
267 recursively_intersect( A[1], tmid, t1, deptha,
268 b, u0, u1, depthb,
269 parameters );
270 }
271 }
272 else
273 if( depthb > 0 )
274 {
275 OldBezier B[2];
276 b.split(0.5, B[0], B[1]);
277 double umid = (u0 + u1)*0.5;
278 depthb--;
279 if( intersect_BB( a, B[0] ) )
280 recursively_intersect( a, t0, t1, deptha,
281 B[0], u0, umid, depthb,
282 parameters );
283 if( intersect_BB( a, B[1] ) )
284 recursively_intersect( a, t0, t1, deptha,
285 B[0], umid, u1, depthb,
286 parameters );
287 }
288 else // Both segments are fully subdivided; now do line segments
289 {
290 double xlk = a.p.back()[X] - a.p[0][X];
291 double ylk = a.p.back()[Y] - a.p[0][Y];
292 double xnm = b.p.back()[X] - b.p[0][X];
293 double ynm = b.p.back()[Y] - b.p[0][Y];
294 double xmk = b.p[0][X] - a.p[0][X];
295 double ymk = b.p[0][Y] - a.p[0][Y];
296 double det = xnm * ylk - ynm * xlk;
297 if( 1.0 + det == 1.0 )
298 return;
299 else
300 {
301 double detinv = 1.0 / det;
302 double s = ( xnm * ymk - ynm *xmk ) * detinv;
303 double t = ( xlk * ymk - ylk * xmk ) * detinv;
304 if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) )
305 return;
306 parameters.push_back(std::pair<double, double>(t0 + s * ( t1 - t0 ),
307 u0 + t * ( u1 - u0 )));
308 }
309 }
310 }
312 inline double log4( double x ) { return log(x)/log(4.); }
314 /*
315 * Wang's theorem is used to estimate the level of subdivision required,
316 * but only if the bounding boxes interfere at the top level.
317 * Assuming there is a possible intersection, recursively_intersect is
318 * used to find all the parameters corresponding to intersection points.
319 * these are then sorted and returned in an array.
320 */
322 double Lmax(Point p) {
323 return std::max(fabs(p[X]), fabs(p[Y]));
324 }
326 unsigned wangs_theorem(OldBezier a) {
327 return 12; // seems a good approximation!
328 double la1 = Lmax( ( a.p[2] - a.p[1] ) - (a.p[1] - a.p[0]) );
329 double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
330 double l0 = std::max(la1, la2);
331 unsigned ra;
332 if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
333 ra = 0;
334 else
335 ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
336 std::cout << ra << std::endl;
337 return ra;
338 }
340 std::vector<std::pair<double, double> > find_intersections( OldBezier a, OldBezier b)
341 {
342 std::vector<std::pair<double, double> > parameters;
343 if( intersect_BB( a, b ) )
344 {
345 recursively_intersect( a, 0., 1., wangs_theorem(a),
346 b, 0., 1., wangs_theorem(b),
347 parameters);
348 }
349 std::sort(parameters.begin(), parameters.end());
350 return parameters;
351 }
352 };
354 /*
355 Local Variables:
356 mode:c++
357 c-file-style:"stroustrup"
358 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
359 indent-tabs-mode:nil
360 fill-column:99
361 End:
362 */
363 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :