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