1 /**************************************************************
2 *
3 * Licensed to the Apache Software Foundation (ASF) under one
4 * or more contributor license agreements. See the NOTICE file
5 * distributed with this work for additional information
6 * regarding copyright ownership. The ASF licenses this file
7 * to you under the Apache License, Version 2.0 (the
8 * "License"); you may not use this file except in compliance
9 * with the License. You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing,
14 * software distributed under the License is distributed on an
15 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16 * KIND, either express or implied. See the License for the
17 * specific language governing permissions and limitations
18 * under the License.
19 *
20 *************************************************************/
21
22
23
24 // MARKER(update_precomp.py): autogen include statement, do not remove
25 #include "precompiled_basegfx.hxx"
26
27 #include <algorithm>
28 #include <iterator>
29 #include <vector>
30 #include <utility>
31
32 #include <math.h>
33
34 #include "bezierclip.hxx"
35 #include "gauss.hxx"
36
37
38
39 // what to test
40 #define WITH_ASSERTIONS
41 //#define WITH_CONVEXHULL_TEST
42 //#define WITH_MULTISUBDIVIDE_TEST
43 //#define WITH_FATLINE_TEST
44 //#define WITH_CALCFOCUS_TEST
45 //#define WITH_SAFEPARAMBASE_TEST
46 //#define WITH_SAFEPARAMS_TEST
47 //#define WITH_SAFEPARAM_DETAILED_TEST
48 //#define WITH_SAFEFOCUSPARAM_CALCFOCUS
49 //#define WITH_SAFEFOCUSPARAM_TEST
50 //#define WITH_SAFEFOCUSPARAM_DETAILED_TEST
51 #define WITH_BEZIERCLIP_TEST
52
53
54
55 // -----------------------------------------------------------------------------
56
57 /* Implementation of the so-called 'Fat-Line Bezier Clipping Algorithm' by Sederberg et al.
58 *
59 * Actual reference is: T. W. Sederberg and T Nishita: Curve
60 * intersection using Bezier clipping. In Computer Aided Design, 22
61 * (9), 1990, pp. 538--549
62 */
63
64 // -----------------------------------------------------------------------------
65
66 /* Misc helper
67 * ===========
68 */
fallFac(int n,int k)69 int fallFac( int n, int k )
70 {
71 #ifdef WITH_ASSERTIONS
72 assert(n>=k); // "For factorials, n must be greater or equal k"
73 assert(n>=0); // "For factorials, n must be positive"
74 assert(k>=0); // "For factorials, k must be positive"
75 #endif
76
77 int res( 1 );
78
79 while( k-- && n ) res *= n--;
80
81 return res;
82 }
83
84 // -----------------------------------------------------------------------------
85
fac(int n)86 int fac( int n )
87 {
88 return fallFac(n, n);
89 }
90
91 // -----------------------------------------------------------------------------
92
93 /* Bezier fat line clipping part
94 * =============================
95 */
96
97 // -----------------------------------------------------------------------------
98
Impl_calcFatLine(FatLine & line,const Bezier & c)99 void Impl_calcFatLine( FatLine& line, const Bezier& c )
100 {
101 // Prepare normalized implicit line
102 // ================================
103
104 // calculate vector orthogonal to p1-p4:
105 line.a = -(c.p0.y - c.p3.y);
106 line.b = (c.p0.x - c.p3.x);
107
108 // normalize
109 const double len( sqrt( line.a*line.a + line.b*line.b ) );
110 if( !tolZero(len) )
111 {
112 line.a /= len;
113 line.b /= len;
114 }
115
116 line.c = -(line.a*c.p0.x + line.b*c.p0.y);
117
118
119 // Determine bounding fat line from it
120 // ===================================
121
122 // calc control point distances
123 const double dP2( calcLineDistance(line.a, line.b, line.c, c.p1.x, c.p1.y ) );
124 const double dP3( calcLineDistance(line.a, line.b, line.c, c.p2.x, c.p2.y ) );
125
126 // calc approximate bounding lines to curve (tight bounds are
127 // possible here, but more expensive to calculate and thus not
128 // worth the overhead)
129 if( dP2 * dP3 > 0.0 )
130 {
131 line.dMin = 3.0/4.0 * ::std::min(0.0, ::std::min(dP2, dP3));
132 line.dMax = 3.0/4.0 * ::std::max(0.0, ::std::max(dP2, dP3));
133 }
134 else
135 {
136 line.dMin = 4.0/9.0 * ::std::min(0.0, ::std::min(dP2, dP3));
137 line.dMax = 4.0/9.0 * ::std::max(0.0, ::std::max(dP2, dP3));
138 }
139 }
140
Impl_calcBounds(Point2D & leftTop,Point2D & rightBottom,const Bezier & c1)141 void Impl_calcBounds( Point2D& leftTop,
142 Point2D& rightBottom,
143 const Bezier& c1 )
144 {
145 leftTop.x = ::std::min( c1.p0.x, ::std::min( c1.p1.x, ::std::min( c1.p2.x, c1.p3.x ) ) );
146 leftTop.y = ::std::min( c1.p0.y, ::std::min( c1.p1.y, ::std::min( c1.p2.y, c1.p3.y ) ) );
147 rightBottom.x = ::std::max( c1.p0.x, ::std::max( c1.p1.x, ::std::max( c1.p2.x, c1.p3.x ) ) );
148 rightBottom.y = ::std::max( c1.p0.y, ::std::max( c1.p1.y, ::std::max( c1.p2.y, c1.p3.y ) ) );
149 }
150
Impl_doBBoxIntersect(const Bezier & c1,const Bezier & c2)151 bool Impl_doBBoxIntersect( const Bezier& c1,
152 const Bezier& c2 )
153 {
154 // calc rectangular boxes from c1 and c2
155 Point2D lt1;
156 Point2D rb1;
157 Point2D lt2;
158 Point2D rb2;
159
160 Impl_calcBounds( lt1, rb1, c1 );
161 Impl_calcBounds( lt2, rb2, c2 );
162
163 if( ::std::min(rb1.x, rb2.x) < ::std::max(lt1.x, lt2.x) ||
164 ::std::min(rb1.y, rb2.y) < ::std::max(lt1.y, lt2.y) )
165 {
166 return false;
167 }
168 else
169 {
170 return true;
171 }
172 }
173
174 /* calculates two t's for the given bernstein control polygon: the first is
175 * the intersection of the min value line with the convex hull from
176 * the left, the second is the intersection of the max value line with
177 * the convex hull from the right.
178 */
Impl_calcSafeParams(double & t1,double & t2,const Polygon2D & rPoly,double lowerYBound,double upperYBound)179 bool Impl_calcSafeParams( double& t1,
180 double& t2,
181 const Polygon2D& rPoly,
182 double lowerYBound,
183 double upperYBound )
184 {
185 // need the convex hull of the control polygon, as this is
186 // guaranteed to completely bound the curve
187 Polygon2D convHull( convexHull(rPoly) );
188
189 // init min and max buffers
190 t1 = 0.0 ;
191 double currLowerT( 1.0 );
192
193 t2 = 1.0;
194 double currHigherT( 0.0 );
195
196 if( convHull.size() <= 1 )
197 return false; // only one point? Then we're done with clipping
198
199 /* now, clip against lower and higher bounds */
200 Point2D p0;
201 Point2D p1;
202
203 bool bIntersection( false );
204
205 for( Polygon2D::size_type i=0; i<convHull.size(); ++i )
206 {
207 // have to check against convHull.size() segments, as the
208 // convex hull is, by definition, closed. Thus, for the
209 // last point, we take the first point as partner.
210 if( i+1 == convHull.size() )
211 {
212 // close the polygon
213 p0 = convHull[i];
214 p1 = convHull[0];
215 }
216 else
217 {
218 p0 = convHull[i];
219 p1 = convHull[i+1];
220 }
221
222 // is the segment in question within or crossing the
223 // horizontal band spanned by lowerYBound and upperYBound? If
224 // not, we've got no intersection. If yes, we maybe don't have
225 // an intersection, but we've got to update the permissible
226 // range, nevertheless. This is because inside lying segments
227 // leads to full range forbidden.
228 if( (tolLessEqual(p0.y, upperYBound) || tolLessEqual(p1.y, upperYBound)) &&
229 (tolGreaterEqual(p0.y, lowerYBound) || tolGreaterEqual(p1.y, lowerYBound)) )
230 {
231 // calc intersection of convex hull segment with
232 // one of the horizontal bounds lines
233 const double r_x( p1.x - p0.x );
234 const double r_y( p1.y - p0.y );
235
236 if( tolZero(r_y) )
237 {
238 // r_y is virtually zero, thus we've got a horizontal
239 // line. Now check whether we maybe coincide with lower or
240 // upper horizonal bound line.
241 if( tolEqual(p0.y, lowerYBound) ||
242 tolEqual(p0.y, upperYBound) )
243 {
244 // yes, simulate intersection then
245 currLowerT = ::std::min(currLowerT, ::std::min(p0.x, p1.x));
246 currHigherT = ::std::max(currHigherT, ::std::max(p0.x, p1.x));
247 }
248 }
249 else
250 {
251 // check against lower and higher bounds
252 // =====================================
253
254 // calc intersection with horizontal dMin line
255 const double currTLow( (lowerYBound - p0.y) * r_x / r_y + p0.x );
256
257 // calc intersection with horizontal dMax line
258 const double currTHigh( (upperYBound - p0.y) * r_x / r_y + p0.x );
259
260 currLowerT = ::std::min(currLowerT, ::std::min(currTLow, currTHigh));
261 currHigherT = ::std::max(currHigherT, ::std::max(currTLow, currTHigh));
262 }
263
264 // set flag that at least one segment is contained or
265 // intersects given horizontal band.
266 bIntersection = true;
267 }
268 }
269
270 #ifndef WITH_SAFEPARAMBASE_TEST
271 // limit intersections found to permissible t parameter range
272 t1 = ::std::max(0.0, currLowerT);
273 t2 = ::std::min(1.0, currHigherT);
274 #endif
275
276 return bIntersection;
277 }
278
279
280 /* calculates two t's for the given bernstein polynomial: the first is
281 * the intersection of the min value line with the convex hull from
282 * the left, the second is the intersection of the max value line with
283 * the convex hull from the right.
284 *
285 * The polynomial coefficients c0 to c3 given to this method
286 * must correspond to t values of 0, 1/3, 2/3 and 1, respectively.
287 */
Impl_calcSafeParams_clip(double & t1,double & t2,const FatLine & bounds,double c0,double c1,double c2,double c3)288 bool Impl_calcSafeParams_clip( double& t1,
289 double& t2,
290 const FatLine& bounds,
291 double c0,
292 double c1,
293 double c2,
294 double c3 )
295 {
296 /* first of all, determine convex hull of c0-c3 */
297 Polygon2D poly(4);
298 poly[0] = Point2D(0, c0);
299 poly[1] = Point2D(1.0/3.0, c1);
300 poly[2] = Point2D(2.0/3.0, c2);
301 poly[3] = Point2D(1, c3);
302
303 #ifndef WITH_SAFEPARAM_DETAILED_TEST
304
305 return Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax );
306
307 #else
308 bool bRet( Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax ) );
309
310 Polygon2D convHull( convexHull( poly ) );
311
312 cout << "# convex hull testing" << endl
313 << "plot [t=0:1] ";
314 cout << " bez("
315 << poly[0].x << ","
316 << poly[1].x << ","
317 << poly[2].x << ","
318 << poly[3].x << ",t),bez("
319 << poly[0].y << ","
320 << poly[1].y << ","
321 << poly[2].y << ","
322 << poly[3].y << ",t), "
323 << "t, " << bounds.dMin << ", "
324 << "t, " << bounds.dMax << ", "
325 << t1 << ", t, "
326 << t2 << ", t, "
327 << "'-' using ($1):($2) title \"control polygon\" with lp, "
328 << "'-' using ($1):($2) title \"convex hull\" with lp" << endl;
329
330 unsigned int k;
331 for( k=0; k<poly.size(); ++k )
332 {
333 cout << poly[k].x << " " << poly[k].y << endl;
334 }
335 cout << poly[0].x << " " << poly[0].y << endl;
336 cout << "e" << endl;
337
338 for( k=0; k<convHull.size(); ++k )
339 {
340 cout << convHull[k].x << " " << convHull[k].y << endl;
341 }
342 cout << convHull[0].x << " " << convHull[0].y << endl;
343 cout << "e" << endl;
344
345 return bRet;
346 #endif
347 }
348
349 // -----------------------------------------------------------------------------
350
Impl_deCasteljauAt(Bezier & part1,Bezier & part2,const Bezier & input,double t)351 void Impl_deCasteljauAt( Bezier& part1,
352 Bezier& part2,
353 const Bezier& input,
354 double t )
355 {
356 // deCasteljau bezier arc, scheme is:
357 //
358 // First row is C_0^n,C_1^n,...,C_n^n
359 // Second row is P_1^n,...,P_n^n
360 // etc.
361 // with P_k^r = (1 - x_s)P_{k-1}^{r-1} + x_s P_k{r-1}
362 //
363 // this results in:
364 //
365 // P1 P2 P3 P4
366 // L1 P2 P3 R4
367 // L2 H R3
368 // L3 R2
369 // L4/R1
370 if( tolZero(t) )
371 {
372 // t is zero -> part2 is input curve, part1 is empty (input.p0, that is)
373 part1.p0.x = part1.p1.x = part1.p2.x = part1.p3.x = input.p0.x;
374 part1.p0.y = part1.p1.y = part1.p2.y = part1.p3.y = input.p0.y;
375 part2 = input;
376 }
377 else if( tolEqual(t, 1.0) )
378 {
379 // t is one -> part1 is input curve, part2 is empty (input.p3, that is)
380 part1 = input;
381 part2.p0.x = part2.p1.x = part2.p2.x = part2.p3.x = input.p3.x;
382 part2.p0.y = part2.p1.y = part2.p2.y = part2.p3.y = input.p3.y;
383 }
384 else
385 {
386 part1.p0.x = input.p0.x; part1.p0.y = input.p0.y;
387 part1.p1.x = (1.0 - t)*part1.p0.x + t*input.p1.x; part1.p1.y = (1.0 - t)*part1.p0.y + t*input.p1.y;
388 const double Hx ( (1.0 - t)*input.p1.x + t*input.p2.x ), Hy ( (1.0 - t)*input.p1.y + t*input.p2.y );
389 part1.p2.x = (1.0 - t)*part1.p1.x + t*Hx; part1.p2.y = (1.0 - t)*part1.p1.y + t*Hy;
390 part2.p3.x = input.p3.x; part2.p3.y = input.p3.y;
391 part2.p2.x = (1.0 - t)*input.p2.x + t*input.p3.x; part2.p2.y = (1.0 - t)*input.p2.y + t*input.p3.y;
392 part2.p1.x = (1.0 - t)*Hx + t*part2.p2.x; part2.p1.y = (1.0 - t)*Hy + t*part2.p2.y;
393 part2.p0.x = (1.0 - t)*part1.p2.x + t*part2.p1.x; part2.p0.y = (1.0 - t)*part1.p2.y + t*part2.p1.y;
394 part1.p3.x = part2.p0.x; part1.p3.y = part2.p0.y;
395 }
396 }
397
398 // -----------------------------------------------------------------------------
399
printCurvesWithSafeRange(const Bezier & c1,const Bezier & c2,double t1_c1,double t2_c1,const Bezier & c2_part,const FatLine & bounds_c2)400 void printCurvesWithSafeRange( const Bezier& c1, const Bezier& c2, double t1_c1, double t2_c1,
401 const Bezier& c2_part, const FatLine& bounds_c2 )
402 {
403 static int offset = 0;
404
405 cout << "# safe param range testing" << endl
406 << "plot [t=0.0:1.0] ";
407
408 // clip safe ranges off c1
409 Bezier c1_part1;
410 Bezier c1_part2;
411 Bezier c1_part3;
412
413 // subdivide at t1_c1
414 Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1_c1 );
415 // subdivide at t2_c1
416 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, t2_c1 );
417
418 // output remaining segment (c1_part1)
419
420 cout << "bez("
421 << c1.p0.x+offset << ","
422 << c1.p1.x+offset << ","
423 << c1.p2.x+offset << ","
424 << c1.p3.x+offset << ",t),bez("
425 << c1.p0.y << ","
426 << c1.p1.y << ","
427 << c1.p2.y << ","
428 << c1.p3.y << ",t), bez("
429 << c2.p0.x+offset << ","
430 << c2.p1.x+offset << ","
431 << c2.p2.x+offset << ","
432 << c2.p3.x+offset << ",t),bez("
433 << c2.p0.y << ","
434 << c2.p1.y << ","
435 << c2.p2.y << ","
436 << c2.p3.y << ",t), "
437 #if 1
438 << "bez("
439 << c1_part1.p0.x+offset << ","
440 << c1_part1.p1.x+offset << ","
441 << c1_part1.p2.x+offset << ","
442 << c1_part1.p3.x+offset << ",t),bez("
443 << c1_part1.p0.y << ","
444 << c1_part1.p1.y << ","
445 << c1_part1.p2.y << ","
446 << c1_part1.p3.y << ",t), "
447 #endif
448 #if 1
449 << "bez("
450 << c2_part.p0.x+offset << ","
451 << c2_part.p1.x+offset << ","
452 << c2_part.p2.x+offset << ","
453 << c2_part.p3.x+offset << ",t),bez("
454 << c2_part.p0.y << ","
455 << c2_part.p1.y << ","
456 << c2_part.p2.y << ","
457 << c2_part.p3.y << ",t), "
458 #endif
459 << "linex("
460 << bounds_c2.a << ","
461 << bounds_c2.b << ","
462 << bounds_c2.c << ",t)+" << offset << ", liney("
463 << bounds_c2.a << ","
464 << bounds_c2.b << ","
465 << bounds_c2.c << ",t) title \"fat line (center)\", linex("
466 << bounds_c2.a << ","
467 << bounds_c2.b << ","
468 << bounds_c2.c-bounds_c2.dMin << ",t)+" << offset << ", liney("
469 << bounds_c2.a << ","
470 << bounds_c2.b << ","
471 << bounds_c2.c-bounds_c2.dMin << ",t) title \"fat line (min) \", linex("
472 << bounds_c2.a << ","
473 << bounds_c2.b << ","
474 << bounds_c2.c-bounds_c2.dMax << ",t)+" << offset << ", liney("
475 << bounds_c2.a << ","
476 << bounds_c2.b << ","
477 << bounds_c2.c-bounds_c2.dMax << ",t) title \"fat line (max) \"" << endl;
478
479 offset += 1;
480 }
481
482 // -----------------------------------------------------------------------------
483
printResultWithFinalCurves(const Bezier & c1,const Bezier & c1_part,const Bezier & c2,const Bezier & c2_part,double t1_c1,double t2_c1)484 void printResultWithFinalCurves( const Bezier& c1, const Bezier& c1_part,
485 const Bezier& c2, const Bezier& c2_part,
486 double t1_c1, double t2_c1 )
487 {
488 static int offset = 0;
489
490 cout << "# final result" << endl
491 << "plot [t=0.0:1.0] ";
492
493 cout << "bez("
494 << c1.p0.x+offset << ","
495 << c1.p1.x+offset << ","
496 << c1.p2.x+offset << ","
497 << c1.p3.x+offset << ",t),bez("
498 << c1.p0.y << ","
499 << c1.p1.y << ","
500 << c1.p2.y << ","
501 << c1.p3.y << ",t), bez("
502 << c1_part.p0.x+offset << ","
503 << c1_part.p1.x+offset << ","
504 << c1_part.p2.x+offset << ","
505 << c1_part.p3.x+offset << ",t),bez("
506 << c1_part.p0.y << ","
507 << c1_part.p1.y << ","
508 << c1_part.p2.y << ","
509 << c1_part.p3.y << ",t), "
510 << " pointmarkx(bez("
511 << c1.p0.x+offset << ","
512 << c1.p1.x+offset << ","
513 << c1.p2.x+offset << ","
514 << c1.p3.x+offset << ","
515 << t1_c1 << "),t), "
516 << " pointmarky(bez("
517 << c1.p0.y << ","
518 << c1.p1.y << ","
519 << c1.p2.y << ","
520 << c1.p3.y << ","
521 << t1_c1 << "),t), "
522 << " pointmarkx(bez("
523 << c1.p0.x+offset << ","
524 << c1.p1.x+offset << ","
525 << c1.p2.x+offset << ","
526 << c1.p3.x+offset << ","
527 << t2_c1 << "),t), "
528 << " pointmarky(bez("
529 << c1.p0.y << ","
530 << c1.p1.y << ","
531 << c1.p2.y << ","
532 << c1.p3.y << ","
533 << t2_c1 << "),t), "
534
535 << "bez("
536 << c2.p0.x+offset << ","
537 << c2.p1.x+offset << ","
538 << c2.p2.x+offset << ","
539 << c2.p3.x+offset << ",t),bez("
540 << c2.p0.y << ","
541 << c2.p1.y << ","
542 << c2.p2.y << ","
543 << c2.p3.y << ",t), "
544 << "bez("
545 << c2_part.p0.x+offset << ","
546 << c2_part.p1.x+offset << ","
547 << c2_part.p2.x+offset << ","
548 << c2_part.p3.x+offset << ",t),bez("
549 << c2_part.p0.y << ","
550 << c2_part.p1.y << ","
551 << c2_part.p2.y << ","
552 << c2_part.p3.y << ",t)" << endl;
553
554 offset += 1;
555 }
556
557 // -----------------------------------------------------------------------------
558
559 /** determine parameter ranges [0,t1) and (t2,1] on c1, where c1 is guaranteed to lie outside c2.
560 Returns false, if the two curves don't even intersect.
561
562 @param t1
563 Range [0,t1) on c1 is guaranteed to lie outside c2
564
565 @param t2
566 Range (t2,1] on c1 is guaranteed to lie outside c2
567
568 @param c1_orig
569 Original curve c1
570
571 @param c1_part
572 Subdivided current part of c1
573
574 @param c2_orig
575 Original curve c2
576
577 @param c2_part
578 Subdivided current part of c2
579 */
Impl_calcClipRange(double & t1,double & t2,const Bezier & c1_orig,const Bezier & c1_part,const Bezier & c2_orig,const Bezier & c2_part)580 bool Impl_calcClipRange( double& t1,
581 double& t2,
582 const Bezier& c1_orig,
583 const Bezier& c1_part,
584 const Bezier& c2_orig,
585 const Bezier& c2_part )
586 {
587 // TODO: Maybe also check fat line orthogonal to P0P3, having P0
588 // and P3 as the extremal points
589
590 if( Impl_doBBoxIntersect(c1_part, c2_part) )
591 {
592 // Calculate fat lines around c1
593 FatLine bounds_c2;
594
595 // must use the subdivided version of c2, since the fat line
596 // algorithm works implicitely with the convex hull bounding
597 // box.
598 Impl_calcFatLine(bounds_c2, c2_part);
599
600 // determine clip positions on c2. Can use original c1 (which
601 // is necessary anyway, to get the t's on the original curve),
602 // since the distance calculations work directly in the
603 // Bernstein polynom parameter domain.
604 if( Impl_calcSafeParams_clip( t1, t2, bounds_c2,
605 calcLineDistance( bounds_c2.a,
606 bounds_c2.b,
607 bounds_c2.c,
608 c1_orig.p0.x,
609 c1_orig.p0.y ),
610 calcLineDistance( bounds_c2.a,
611 bounds_c2.b,
612 bounds_c2.c,
613 c1_orig.p1.x,
614 c1_orig.p1.y ),
615 calcLineDistance( bounds_c2.a,
616 bounds_c2.b,
617 bounds_c2.c,
618 c1_orig.p2.x,
619 c1_orig.p2.y ),
620 calcLineDistance( bounds_c2.a,
621 bounds_c2.b,
622 bounds_c2.c,
623 c1_orig.p3.x,
624 c1_orig.p3.y ) ) )
625 {
626 //printCurvesWithSafeRange(c1_orig, c2_orig, t1, t2, c2_part, bounds_c2);
627
628 // they do intersect
629 return true;
630 }
631 }
632
633 // they don't intersect: nothing to do
634 return false;
635 }
636
637 // -----------------------------------------------------------------------------
638
639 /* Tangent intersection part
640 * =========================
641 */
642
643 // -----------------------------------------------------------------------------
644
Impl_calcFocus(Bezier & res,const Bezier & c)645 void Impl_calcFocus( Bezier& res, const Bezier& c )
646 {
647 // arbitrary small value, for now
648 // TODO: find meaningful value
649 const double minPivotValue( 1.0e-20 );
650
651 Point2D::value_type fMatrix[6];
652 Point2D::value_type fRes[2];
653
654 // calc new curve from hodograph, c and linear blend
655
656 // Coefficients for derivative of c are (C_i=n(C_{i+1} - C_i)):
657 //
658 // 3(P1 - P0), 3(P2 - P1), 3(P3 - P2) (bezier curve of degree 2)
659 //
660 // The hodograph is then (bezier curve of 2nd degree is P0(1-t)^2 + 2P1(1-t)t + P2t^2):
661 //
662 // 3(P1 - P0)(1-t)^2 + 6(P2 - P1)(1-t)t + 3(P3 - P2)t^2
663 //
664 // rotate by 90 degrees: x=-y, y=x and you get the normal vector function N(t):
665 //
666 // x(t) = -(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2)
667 // y(t) = 3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2
668 //
669 // Now, the focus curve is defined to be F(t)=P(t) + c(t)N(t),
670 // where P(t) is the original curve, and c(t)=c0(1-t) + c1 t
671 //
672 // This results in the following expression for F(t):
673 //
674 // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
675 // (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2)
676 //
677 // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1.t)t^2 + P3.y t^3 +
678 // (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2)
679 //
680 // As a heuristic, we set F(0)=F(1) (thus, the curve is closed and _tends_ to be small):
681 //
682 // For F(0), the following results:
683 //
684 // x(0) = P0.x - c0 3(P1.y - P0.y)
685 // y(0) = P0.y + c0 3(P1.x - P0.x)
686 //
687 // For F(1), the following results:
688 //
689 // x(1) = P3.x - c1 3(P3.y - P2.y)
690 // y(1) = P3.y + c1 3(P3.x - P2.x)
691 //
692 // Reorder, collect and substitute into F(0)=F(1):
693 //
694 // P0.x - c0 3(P1.y - P0.y) = P3.x - c1 3(P3.y - P2.y)
695 // P0.y + c0 3(P1.x - P0.x) = P3.y + c1 3(P3.x - P2.x)
696 //
697 // which yields
698 //
699 // (P0.y - P1.y)c0 + (P3.y - P2.y)c1 = (P3.x - P0.x)/3
700 // (P1.x - P0.x)c0 + (P2.x - P3.x)c1 = (P3.y - P0.y)/3
701 //
702
703 // so, this is what we calculate here (determine c0 and c1):
704 fMatrix[0] = c.p1.x - c.p0.x;
705 fMatrix[1] = c.p2.x - c.p3.x;
706 fMatrix[2] = (c.p3.y - c.p0.y)/3.0;
707 fMatrix[3] = c.p0.y - c.p1.y;
708 fMatrix[4] = c.p3.y - c.p2.y;
709 fMatrix[5] = (c.p3.x - c.p0.x)/3.0;
710
711 // TODO: determine meaningful value for
712 if( !solve(fMatrix, 2, 3, fRes, minPivotValue) )
713 {
714 // TODO: generate meaningful values here
715 // singular or nearly singular system -- use arbitrary
716 // values for res
717 fRes[0] = 0.0;
718 fRes[1] = 1.0;
719
720 cerr << "Matrix singular!" << endl;
721 }
722
723 // now, the reordered and per-coefficient collected focus curve is
724 // the following third degree bezier curve F(t):
725 //
726 // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
727 // (c0(1-t) + c1 t)(3(P1.y - P0.y)(1-t)^2 + 6(P2.y - P1.y)(1-t)t + 3(P3.y - P2.y)t^2)
728 // = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
729 // (3c0P1.y(1-t)^3 - 3c0P0.y(1-t)^3 + 6c0P2.y(1-t)^2t - 6c0P1.y(1-t)^2t +
730 // 3c0P3.y(1-t)t^2 - 3c0P2.y(1-t)t^2 +
731 // 3c1P1.y(1-t)^2t - 3c1P0.y(1-t)^2t + 6c1P2.y(1-t)t^2 - 6c1P1.y(1-t)t^2 +
732 // 3c1P3.yt^3 - 3c1P2.yt^3)
733 // = (P0.x - 3 c0 P1.y + 3 c0 P0.y)(1-t)^3 +
734 // 3(P1.x - c1 P1.y + c1 P0.y - 2 c0 P2.y + 2 c0 P1.y)(1-t)^2t +
735 // 3(P2.x - 2 c1 P2.y + 2 c1 P1.y - c0 P3.y + c0 P2.y)(1-t)t^2 +
736 // (P3.x - 3 c1 P3.y + 3 c1 P2.y)t^3
737 // = (P0.x - 3 c0(P1.y - P0.y))(1-t)^3 +
738 // 3(P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y))(1-t)^2t +
739 // 3(P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y))(1-t)t^2 +
740 // (P3.x - 3 c1(P3.y - P2.y))t^3
741 //
742 // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 +
743 // (c0(1-t) + c1 t)(3(P1.x - P0.x)(1-t)^2 + 6(P2.x - P1.x)(1-t)t + 3(P3.x - P2.x)t^2)
744 // = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 +
745 // 3c0(P1.x - P0.x)(1-t)^3 + 6c0(P2.x - P1.x)(1-t)^2t + 3c0(P3.x - P2.x)(1-t)t^2 +
746 // 3c1(P1.x - P0.x)(1-t)^2t + 6c1(P2.x - P1.x)(1-t)t^2 + 3c1(P3.x - P2.x)t^3
747 // = (P0.y + 3 c0 (P1.x - P0.x))(1-t)^3 +
748 // 3(P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x))(1-t)^2t +
749 // 3(P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x))(1-t)t^2 +
750 // (P3.y + 3 c1 (P3.x - P2.x))t^3
751 //
752 // Therefore, the coefficients F0 to F3 of the focus curve are:
753 //
754 // F0.x = (P0.x - 3 c0(P1.y - P0.y)) F0.y = (P0.y + 3 c0 (P1.x - P0.x))
755 // F1.x = (P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y)) F1.y = (P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x))
756 // F2.x = (P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y)) F2.y = (P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x))
757 // F3.x = (P3.x - 3 c1(P3.y - P2.y)) F3.y = (P3.y + 3 c1 (P3.x - P2.x))
758 //
759 res.p0.x = c.p0.x - 3*fRes[0]*(c.p1.y - c.p0.y);
760 res.p1.x = c.p1.x - fRes[1]*(c.p1.y - c.p0.y) - 2*fRes[0]*(c.p2.y - c.p1.y);
761 res.p2.x = c.p2.x - 2*fRes[1]*(c.p2.y - c.p1.y) - fRes[0]*(c.p3.y - c.p2.y);
762 res.p3.x = c.p3.x - 3*fRes[1]*(c.p3.y - c.p2.y);
763
764 res.p0.y = c.p0.y + 3*fRes[0]*(c.p1.x - c.p0.x);
765 res.p1.y = c.p1.y + 2*fRes[0]*(c.p2.x - c.p1.x) + fRes[1]*(c.p1.x - c.p0.x);
766 res.p2.y = c.p2.y + fRes[0]*(c.p3.x - c.p2.x) + 2*fRes[1]*(c.p2.x - c.p1.x);
767 res.p3.y = c.p3.y + 3*fRes[1]*(c.p3.x - c.p2.x);
768 }
769
770 // -----------------------------------------------------------------------------
771
Impl_calcSafeParams_focus(double & t1,double & t2,const Bezier & curve,const Bezier & focus)772 bool Impl_calcSafeParams_focus( double& t1,
773 double& t2,
774 const Bezier& curve,
775 const Bezier& focus )
776 {
777 // now, we want to determine which normals of the original curve
778 // P(t) intersect with the focus curve F(t). The condition for
779 // this statement is P'(t)(P(t) - F) = 0, i.e. hodograph P'(t) and
780 // line through P(t) and F are perpendicular.
781 // If you expand this equation, you end up with something like
782 //
783 // (\sum_{i=0}^n (P_i - F)B_i^n(t))^T (\sum_{j=0}^{n-1} n(P_{j+1} - P_j)B_j^{n-1}(t))
784 //
785 // Multiplying that out (as the scalar product is linear, we can
786 // extract some terms) yields:
787 //
788 // (P_i - F)^T n(P_{j+1} - P_j) B_i^n(t)B_j^{n-1}(t) + ...
789 //
790 // If we combine the B_i^n(t)B_j^{n-1}(t) product, we arrive at a
791 // Bernstein polynomial of degree 2n-1, as
792 //
793 // \binom{n}{i}(1-t)^{n-i}t^i) \binom{n-1}{j}(1-t)^{n-1-j}t^j) =
794 // \binom{n}{i}\binom{n-1}{j}(1-t)^{2n-1-i-j}t^{i+j}
795 //
796 // Thus, with the defining equation for a 2n-1 degree Bernstein
797 // polynomial
798 //
799 // \sum_{i=0}^{2n-1} d_i B_i^{2n-1}(t)
800 //
801 // the d_i are calculated as follows:
802 //
803 // d_i = \sum_{j+k=i, j\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{j}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_j - F)
804 //
805 //
806 // Okay, but F is now not a single point, but itself a curve
807 // F(u). Thus, for every value of u, we get a different 2n-1
808 // bezier curve from the above equation. Therefore, we have a
809 // tensor product bezier patch, with the following defining
810 // equation:
811 //
812 // d(t,u) = \sum_{i=0}^{2n-1} \sum_{j=0}^m B_i^{2n-1}(t) B_j^{m}(u) d_{ij}, where
813 // d_{ij} = \sum_{k+l=i, l\in\{0,...,n\}, k\in\{0,...,n-1\}} \frac{\binom{n}{l}\binom{n-1}{k}}{\binom{2n-1}{i}} n (P_{k+1} - P_k)^T(P_l - F_j)
814 //
815 // as above, only that now F is one of the focus' control points.
816 //
817 // Note the difference in the binomial coefficients to the
818 // reference paper, these formulas most probably contained a typo.
819 //
820 // To determine, where D(t,u) is _not_ zero (these are the parts
821 // of the curve that don't share normals with the focus and can
822 // thus be safely clipped away), we project D(u,t) onto the
823 // (d(t,u), t) plane, determine the convex hull there and proceed
824 // as for the curve intersection part (projection is orthogonal to
825 // u axis, thus simply throw away u coordinate).
826 //
827 // \fallfac are so-called falling factorials (see Concrete
828 // Mathematics, p. 47 for a definition).
829 //
830
831 // now, for tensor product bezier curves, the convex hull property
832 // holds, too. Thus, we simply project the control points (t_{ij},
833 // u_{ij}, d_{ij}) onto the (t,d) plane and calculate the
834 // intersections of the convex hull with the t axis, as for the
835 // bezier clipping case.
836
837 //
838 // calc polygon of control points (t_{ij}, d_{ij}):
839 //
840 const int n( 3 ); // cubic bezier curves, as a matter of fact
841 const int i_card( 2*n );
842 const int j_card( n + 1 );
843 const int k_max( n-1 );
844 Polygon2D controlPolygon( i_card*j_card ); // vector of (t_{ij}, d_{ij}) in row-major order
845
846 int i, j, k, l; // variable notation from formulas above and Sederberg article
847 Point2D::value_type d;
848 for( i=0; i<i_card; ++i )
849 {
850 for( j=0; j<j_card; ++j )
851 {
852 // calc single d_{ij} sum:
853 for( d=0.0, k=::std::max(0,i-n); k<=k_max && k<=i; ++k )
854 {
855 l = i - k; // invariant: k + l = i
856 assert(k>=0 && k<=n-1); // k \in {0,...,n-1}
857 assert(l>=0 && l<=n); // l \in {0,...,n}
858
859 // TODO: find, document and assert proper limits for n and int's max_val.
860 // This becomes important should anybody wants to use
861 // this code for higher-than-cubic beziers
862 d += static_cast<double>(fallFac(n,l)*fallFac(n-1,k)*fac(i)) /
863 static_cast<double>(fac(l)*fac(k) * fallFac(2*n-1,i)) * n *
864 ( (curve[k+1].x - curve[k].x)*(curve[l].x - focus[j].x) + // dot product here
865 (curve[k+1].y - curve[k].y)*(curve[l].y - focus[j].y) );
866 }
867
868 // Note that the t_{ij} values are evenly spaced on the
869 // [0,1] interval, thus t_{ij}=i/(2n-1)
870 controlPolygon[ i*j_card + j ] = Point2D( i/(2.0*n-1.0), d );
871 }
872 }
873
874 #ifndef WITH_SAFEFOCUSPARAM_DETAILED_TEST
875
876 // calc safe parameter range, to determine [0,t1] and [t2,1] where
877 // no zero crossing is guaranteed.
878 return Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 );
879
880 #else
881 bool bRet( Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 ) );
882
883 Polygon2D convHull( convexHull( controlPolygon ) );
884
885 cout << "# convex hull testing (focus)" << endl
886 << "plot [t=0:1] ";
887 cout << "'-' using ($1):($2) title \"control polygon\" with lp, "
888 << "'-' using ($1):($2) title \"convex hull\" with lp" << endl;
889
890 unsigned int count;
891 for( count=0; count<controlPolygon.size(); ++count )
892 {
893 cout << controlPolygon[count].x << " " << controlPolygon[count].y << endl;
894 }
895 cout << controlPolygon[0].x << " " << controlPolygon[0].y << endl;
896 cout << "e" << endl;
897
898 for( count=0; count<convHull.size(); ++count )
899 {
900 cout << convHull[count].x << " " << convHull[count].y << endl;
901 }
902 cout << convHull[0].x << " " << convHull[0].y << endl;
903 cout << "e" << endl;
904
905 return bRet;
906 #endif
907 }
908
909 // -----------------------------------------------------------------------------
910
911 /** Calc all values t_i on c1, for which safeRanges functor does not
912 give a safe range on c1 and c2.
913
914 This method is the workhorse of the bezier clipping. Because c1
915 and c2 must be alternatingly tested against each other (first
916 determine safe parameter interval on c1 with regard to c2, then
917 the other way around), we call this method recursively with c1 and
918 c2 swapped.
919
920 @param result
921 Output iterator where the final t values are added to. If curves
922 don't intersect, nothing is added.
923
924 @param delta
925 Maximal allowed distance to true critical point (measured in the
926 original curve's coordinate system)
927
928 @param safeRangeFunctor
929 Functor object, that must provide the following operator():
930 bool safeRangeFunctor( double& t1,
931 double& t2,
932 const Bezier& c1_orig,
933 const Bezier& c1_part,
934 const Bezier& c2_orig,
935 const Bezier& c2_part );
936 This functor must calculate the safe ranges [0,t1] and [t2,1] on
937 c1_orig, where c1_orig is 'safe' from c2_part. If the whole
938 c1_orig is safe, false must be returned, true otherwise.
939 */
Impl_applySafeRanges_rec(::std::back_insert_iterator<::std::vector<::std::pair<double,double>>> & result,double delta,const Functor & safeRangeFunctor,int recursionLevel,const Bezier & c1_orig,const Bezier & c1_part,double last_t1_c1,double last_t2_c1,const Bezier & c2_orig,const Bezier & c2_part,double last_t1_c2,double last_t2_c2)940 template <class Functor> void Impl_applySafeRanges_rec( ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > >& result,
941 double delta,
942 const Functor& safeRangeFunctor,
943 int recursionLevel,
944 const Bezier& c1_orig,
945 const Bezier& c1_part,
946 double last_t1_c1,
947 double last_t2_c1,
948 const Bezier& c2_orig,
949 const Bezier& c2_part,
950 double last_t1_c2,
951 double last_t2_c2 )
952 {
953 // check end condition
954 // ===================
955
956 // TODO: tidy up recursion handling. maybe put everything in a
957 // struct and swap that here at method entry
958
959 // TODO: Implement limit on recursion depth. Should that limit be
960 // reached, chances are that we're on a higher-order tangency. For
961 // this case, AW proposed to take the middle of the current
962 // interval, and to correct both curve's tangents at that new
963 // endpoint to be equal. That virtually generates a first-order
964 // tangency, and justifies to return a single intersection
965 // point. Otherwise, inside/outside test might fail here.
966
967 for( int i=0; i<recursionLevel; ++i ) cerr << " ";
968 if( recursionLevel % 2 )
969 {
970 cerr << "level: " << recursionLevel
971 << " t: "
972 << last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0
973 << ", c1: " << last_t1_c2 << " " << last_t2_c2
974 << ", c2: " << last_t1_c1 << " " << last_t2_c1
975 << endl;
976 }
977 else
978 {
979 cerr << "level: " << recursionLevel
980 << " t: "
981 << last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0
982 << ", c1: " << last_t1_c1 << " " << last_t2_c1
983 << ", c2: " << last_t1_c2 << " " << last_t2_c2
984 << endl;
985 }
986
987 // refine solution
988 // ===============
989
990 double t1_c1, t2_c1;
991
992 // Note: we first perform the clipping and only test for precision
993 // sufficiency afterwards, since we want to exploit the fact that
994 // Impl_calcClipRange returns false if the curves don't
995 // intersect. We would have to check that separately for the end
996 // condition, otherwise.
997
998 // determine safe range on c1_orig
999 if( safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part ) )
1000 {
1001 // now, t1 and t2 are calculated on the original curve
1002 // (but against a fat line calculated from the subdivided
1003 // c2, namely c2_part). If the [t1,t2] range is outside
1004 // our current [last_t1,last_t2] range, we're done in this
1005 // branch - the curves no longer intersect.
1006 if( tolLessEqual(t1_c1, last_t2_c1) && tolGreaterEqual(t2_c1, last_t1_c1) )
1007 {
1008 // As noted above, t1 and t2 are calculated on the
1009 // original curve, but against a fat line
1010 // calculated from the subdivided c2, namely
1011 // c2_part. Our domain to work on is
1012 // [last_t1,last_t2], on the other hand, so values
1013 // of [t1,t2] outside that range are irrelevant
1014 // here. Clip range appropriately.
1015 t1_c1 = ::std::max(t1_c1, last_t1_c1);
1016 t2_c1 = ::std::min(t2_c1, last_t2_c1);
1017
1018 // TODO: respect delta
1019 // for now, end condition is just a fixed threshold on the t's
1020
1021 // check end condition
1022 // ===================
1023
1024 #if 1
1025 if( fabs(last_t2_c1 - last_t1_c1) < 0.0001 &&
1026 fabs(last_t2_c2 - last_t1_c2) < 0.0001 )
1027 #else
1028 if( fabs(last_t2_c1 - last_t1_c1) < 0.01 &&
1029 fabs(last_t2_c2 - last_t1_c2) < 0.01 )
1030 #endif
1031 {
1032 // done. Add to result
1033 if( recursionLevel % 2 )
1034 {
1035 // uneven level: have to swap the t's, since curves are swapped, too
1036 *result++ = ::std::make_pair( last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0,
1037 last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0 );
1038 }
1039 else
1040 {
1041 *result++ = ::std::make_pair( last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0,
1042 last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0 );
1043 }
1044
1045 #if 0
1046 //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, last_t1_c1, last_t2_c1 );
1047 printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, t1_c1, t2_c1 );
1048 #else
1049 // calc focus curve of c2
1050 Bezier focus;
1051 Impl_calcFocus(focus, c2_part); // need to use subdivided c2
1052
1053 safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part );
1054
1055 //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, t1_c1, t2_c1 );
1056 printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, last_t1_c1, last_t2_c1 );
1057 #endif
1058 }
1059 else
1060 {
1061 // heuristic: if parameter range is not reduced by at least
1062 // 20%, subdivide longest curve, and clip shortest against
1063 // both parts of longest
1064 // if( (last_t2_c1 - last_t1_c1 - t2_c1 + t1_c1) / (last_t2_c1 - last_t1_c1) < 0.2 )
1065 if( false )
1066 {
1067 // subdivide and descend
1068 // =====================
1069
1070 Bezier part1;
1071 Bezier part2;
1072
1073 double intervalMiddle;
1074
1075 if( last_t2_c1 - last_t1_c1 > last_t2_c2 - last_t1_c2 )
1076 {
1077 // subdivide c1
1078 // ============
1079
1080 intervalMiddle = last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0;
1081
1082 // subdivide at the middle of the interval (as
1083 // we're not subdividing on the original
1084 // curve, this simply amounts to subdivision
1085 // at 0.5)
1086 Impl_deCasteljauAt( part1, part2, c1_part, 0.5 );
1087
1088 // and descend recursively with swapped curves
1089 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1090 c2_orig, c2_part, last_t1_c2, last_t2_c2,
1091 c1_orig, part1, last_t1_c1, intervalMiddle );
1092
1093 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1094 c2_orig, c2_part, last_t1_c2, last_t2_c2,
1095 c1_orig, part2, intervalMiddle, last_t2_c1 );
1096 }
1097 else
1098 {
1099 // subdivide c2
1100 // ============
1101
1102 intervalMiddle = last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0;
1103
1104 // subdivide at the middle of the interval (as
1105 // we're not subdividing on the original
1106 // curve, this simply amounts to subdivision
1107 // at 0.5)
1108 Impl_deCasteljauAt( part1, part2, c2_part, 0.5 );
1109
1110 // and descend recursively with swapped curves
1111 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1112 c2_orig, part1, last_t1_c2, intervalMiddle,
1113 c1_orig, c1_part, last_t1_c1, last_t2_c1 );
1114
1115 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1116 c2_orig, part2, intervalMiddle, last_t2_c2,
1117 c1_orig, c1_part, last_t1_c1, last_t2_c1 );
1118 }
1119 }
1120 else
1121 {
1122 // apply calculated clip
1123 // =====================
1124
1125 // clip safe ranges off c1_orig
1126 Bezier c1_part1;
1127 Bezier c1_part2;
1128 Bezier c1_part3;
1129
1130 // subdivide at t1_c1
1131 Impl_deCasteljauAt( c1_part1, c1_part2, c1_orig, t1_c1 );
1132
1133 // subdivide at t2_c1. As we're working on
1134 // c1_part2 now, we have to adapt t2_c1 since
1135 // we're no longer in the original parameter
1136 // interval. This is based on the following
1137 // assumption: t2_new = (t2-t1)/(1-t1), which
1138 // relates the t2 value into the new parameter
1139 // range [0,1] of c1_part2.
1140 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2_c1-t1_c1)/(1.0-t1_c1) );
1141
1142 // descend with swapped curves and c1_part1 as the
1143 // remaining (middle) segment
1144 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1145 c2_orig, c2_part, last_t1_c2, last_t2_c2,
1146 c1_orig, c1_part1, t1_c1, t2_c1 );
1147 }
1148 }
1149 }
1150 }
1151 }
1152
1153 // -----------------------------------------------------------------------------
1154
1155 struct ClipBezierFunctor
1156 {
operator ()ClipBezierFunctor1157 bool operator()( double& t1_c1,
1158 double& t2_c1,
1159 const Bezier& c1_orig,
1160 const Bezier& c1_part,
1161 const Bezier& c2_orig,
1162 const Bezier& c2_part ) const
1163 {
1164 return Impl_calcClipRange( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part );
1165 }
1166 };
1167
1168 // -----------------------------------------------------------------------------
1169
1170 struct BezierTangencyFunctor
1171 {
operator ()BezierTangencyFunctor1172 bool operator()( double& t1_c1,
1173 double& t2_c1,
1174 const Bezier& c1_orig,
1175 const Bezier& c1_part,
1176 const Bezier& c2_orig,
1177 const Bezier& c2_part ) const
1178 {
1179 // calc focus curve of c2
1180 Bezier focus;
1181 Impl_calcFocus(focus, c2_part); // need to use subdivided c2
1182 // here, as the whole curve is
1183 // used for focus calculation
1184
1185 // determine safe range on c1_orig
1186 bool bRet( Impl_calcSafeParams_focus( t1_c1, t2_c1,
1187 c1_orig, // use orig curve here, need t's on original curve
1188 focus ) );
1189
1190 cerr << "range: " << t2_c1 - t1_c1 << ", ret: " << bRet << endl;
1191
1192 return bRet;
1193 }
1194 };
1195
1196 // -----------------------------------------------------------------------------
1197
1198 /** Perform a bezier clip (curve against curve)
1199
1200 @param result
1201 Output iterator where the final t values are added to. This
1202 iterator will remain empty, if there are no intersections.
1203
1204 @param delta
1205 Maximal allowed distance to true intersection (measured in the
1206 original curve's coordinate system)
1207 */
clipBezier(::std::back_insert_iterator<::std::vector<::std::pair<double,double>>> & result,double delta,const Bezier & c1,const Bezier & c2)1208 void clipBezier( ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > >& result,
1209 double delta,
1210 const Bezier& c1,
1211 const Bezier& c2 )
1212 {
1213 #if 0
1214 // first of all, determine list of collinear normals. Collinear
1215 // normals typically separate two intersections, thus, subdivide
1216 // at all collinear normal's t values beforehand. This will cater
1217 // for tangent intersections, where two or more intersections are
1218 // infinitesimally close together.
1219
1220 // TODO: evaluate effects of higher-than-second-order
1221 // tangencies. Sederberg et al. state that collinear normal
1222 // algorithm then degrades quickly.
1223
1224 ::std::vector< ::std::pair<double,double> > results;
1225 ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > > ii(results);
1226
1227 Impl_calcCollinearNormals( ii, delta, 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
1228
1229 // As Sederberg's collinear normal theorem is only sufficient, not
1230 // necessary for two intersections left and right, we've to test
1231 // all segments generated by the collinear normal algorithm
1232 // against each other. In other words, if the two curves are both
1233 // divided in a left and a right part, the collinear normal
1234 // theorem does _not_ state that the left part of curve 1 does not
1235 // e.g. intersect with the right part of curve 2.
1236
1237 // divide c1 and c2 at collinear normal intersection points
1238 ::std::vector< Bezier > c1_segments( results.size()+1 );
1239 ::std::vector< Bezier > c2_segments( results.size()+1 );
1240 Bezier c1_remainder( c1 );
1241 Bezier c2_remainder( c2 );
1242 unsigned int i;
1243 for( i=0; i<results.size(); ++i )
1244 {
1245 Bezier c1_part2;
1246 Impl_deCasteljauAt( c1_segments[i], c1_part2, c1_remainder, results[i].first );
1247 c1_remainder = c1_part2;
1248
1249 Bezier c2_part2;
1250 Impl_deCasteljauAt( c2_segments[i], c2_part2, c2_remainder, results[i].second );
1251 c2_remainder = c2_part2;
1252 }
1253 c1_segments[i] = c1_remainder;
1254 c2_segments[i] = c2_remainder;
1255
1256 // now, c1/c2_segments contain all segments, then
1257 // clip every resulting segment against every other
1258 unsigned int c1_curr, c2_curr;
1259 for( c1_curr=0; c1_curr<c1_segments.size(); ++c1_curr )
1260 {
1261 for( c2_curr=0; c2_curr<c2_segments.size(); ++c2_curr )
1262 {
1263 if( c1_curr != c2_curr )
1264 {
1265 Impl_clipBezier_rec(result, delta, 0,
1266 c1_segments[c1_curr], c1_segments[c1_curr],
1267 0.0, 1.0,
1268 c2_segments[c2_curr], c2_segments[c2_curr],
1269 0.0, 1.0);
1270 }
1271 }
1272 }
1273 #else
1274 Impl_applySafeRanges_rec( result, delta, BezierTangencyFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
1275 //Impl_applySafeRanges_rec( result, delta, ClipBezierFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
1276 #endif
1277 // that's it, boys'n'girls!
1278 }
1279
main(int argc,const char * argv[])1280 int main(int argc, const char *argv[])
1281 {
1282 double curr_Offset( 0 );
1283 unsigned int i,j,k;
1284
1285 Bezier someCurves[] =
1286 {
1287 // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.0)},
1288 // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)},
1289 // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)}
1290 // {Point2D(0.25+1,0.5),Point2D(0.25+1,0.708333),Point2D(0.423611+1,0.916667),Point2D(0.770833+1,0.980324)},
1291 // {Point2D(0.0+1,0.0),Point2D(0.0+1,1.0),Point2D(1.0+1,1.0),Point2D(1.0+1,0.5)}
1292
1293 // tangency1
1294 // {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)},
1295 // {Point2D(0.0,1.0),Point2D(0.1,1.0),Point2D(0.4,1.0),Point2D(0.5,1.0)}
1296
1297 // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)},
1298 // {Point2D(0.60114,0.933091),Point2D(0.69461,0.969419),Point2D(0.80676,0.992976),Point2D(0.93756,0.998663)}
1299 // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)},
1300 // {Point2D(0.62712,0.828427),Point2D(0.76305,0.828507),Point2D(0.88555,0.77312),Point2D(0.95069,0.67325)}
1301
1302 // clipping1
1303 // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
1304 // {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)}
1305
1306 // tangency2
1307 // {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)},
1308 // {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)}
1309
1310 // tangency3
1311 // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)},
1312 // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)}
1313
1314 // tangency4
1315 // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)},
1316 // {Point2D(0.26,0.4),Point2D(0.25,0.5),Point2D(0.25,0.5),Point2D(0.26,0.6)}
1317
1318 // tangency5
1319 // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
1320 // {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)}
1321
1322 // tangency6
1323 // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
1324 // {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)}
1325
1326 // tangency7
1327 // {Point2D(2.505,0.0),Point2D(2.505+4.915,4.300),Point2D(2.505+3.213,10.019),Point2D(2.505-2.505,10.255)},
1328 // {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)}
1329
1330 // tangency Sederberg example
1331 {Point2D(2.505,0.0),Point2D(2.505+4.915,4.300),Point2D(2.505+3.213,10.019),Point2D(2.505-2.505,10.255)},
1332 {Point2D(5.33+9.311,0.0),Point2D(5.33+9.311-13.279,4.205),Point2D(5.33+9.311-10.681,9.119),Point2D(5.33+9.311-2.603,10.254)}
1333
1334 // clipping2
1335 // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)},
1336 // {Point2D(0.2575,0.4),Point2D(0.2475,0.5),Point2D(0.2475,0.5),Point2D(0.2575,0.6)}
1337
1338 // {Point2D(0.0,0.1),Point2D(0.2,3.5),Point2D(1.0,-2.5),Point2D(1.1,1.2)},
1339 // {Point2D(0.0,1.0),Point2D(3.5,0.9),Point2D(-2.5,0.1),Point2D(1.1,0.2)}
1340 // {Point2D(0.0,0.1),Point2D(0.2,3.0),Point2D(1.0,-2.0),Point2D(1.1,1.2)},
1341 // {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)}
1342 // {Point2D(0.0,1.0),Point2D(3.0,0.9),Point2D(-2.0,0.1),Point2D(1.1,0.2)}
1343 // {Point2D(0.0,4.0),Point2D(0.1,5.0),Point2D(0.9,5.0),Point2D(1.0,4.0)},
1344 // {Point2D(0.0,0.0),Point2D(0.1,0.5),Point2D(0.9,0.5),Point2D(1.0,0.0)},
1345 // {Point2D(0.0,0.1),Point2D(0.1,1.5),Point2D(0.9,1.5),Point2D(1.0,0.1)},
1346 // {Point2D(0.0,-4.0),Point2D(0.1,-5.0),Point2D(0.9,-5.0),Point2D(1.0,-4.0)}
1347 };
1348
1349 // output gnuplot setup
1350 cout << "#!/usr/bin/gnuplot -persist" << endl
1351 << "#" << endl
1352 << "# automatically generated by bezierclip, don't change!" << endl
1353 << "#" << endl
1354 << "set parametric" << endl
1355 << "bez(p,q,r,s,t) = p*(1-t)**3+q*3*(1-t)**2*t+r*3*(1-t)*t**2+s*t**3" << endl
1356 << "bezd(p,q,r,s,t) = 3*(q-p)*(1-t)**2+6*(r-q)*(1-t)*t+3*(s-r)*t**2" << endl
1357 << "pointmarkx(c,t) = c-0.03*t" << endl
1358 << "pointmarky(c,t) = c+0.03*t" << endl
1359 << "linex(a,b,c,t) = a*-c + t*-b" << endl
1360 << "liney(a,b,c,t) = b*-c + t*a" << endl << endl
1361 << "# end of setup" << endl << endl;
1362
1363 #ifdef WITH_CONVEXHULL_TEST
1364 // test convex hull algorithm
1365 const double convHull_xOffset( curr_Offset );
1366 curr_Offset += 20;
1367 cout << "# convex hull testing" << endl
1368 << "plot [t=0:1] ";
1369 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1370 {
1371 Polygon2D aTestPoly(4);
1372 aTestPoly[0] = someCurves[i].p0;
1373 aTestPoly[1] = someCurves[i].p1;
1374 aTestPoly[2] = someCurves[i].p2;
1375 aTestPoly[3] = someCurves[i].p3;
1376
1377 aTestPoly[0].x += convHull_xOffset;
1378 aTestPoly[1].x += convHull_xOffset;
1379 aTestPoly[2].x += convHull_xOffset;
1380 aTestPoly[3].x += convHull_xOffset;
1381
1382 cout << " bez("
1383 << aTestPoly[0].x << ","
1384 << aTestPoly[1].x << ","
1385 << aTestPoly[2].x << ","
1386 << aTestPoly[3].x << ",t),bez("
1387 << aTestPoly[0].y << ","
1388 << aTestPoly[1].y << ","
1389 << aTestPoly[2].y << ","
1390 << aTestPoly[3].y << ",t), '-' using ($1):($2) title \"convex hull " << i << "\" with lp";
1391
1392 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1393 cout << ",\\" << endl;
1394 else
1395 cout << endl;
1396 }
1397 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1398 {
1399 Polygon2D aTestPoly(4);
1400 aTestPoly[0] = someCurves[i].p0;
1401 aTestPoly[1] = someCurves[i].p1;
1402 aTestPoly[2] = someCurves[i].p2;
1403 aTestPoly[3] = someCurves[i].p3;
1404
1405 aTestPoly[0].x += convHull_xOffset;
1406 aTestPoly[1].x += convHull_xOffset;
1407 aTestPoly[2].x += convHull_xOffset;
1408 aTestPoly[3].x += convHull_xOffset;
1409
1410 Polygon2D convHull( convexHull(aTestPoly) );
1411
1412 for( k=0; k<convHull.size(); ++k )
1413 {
1414 cout << convHull[k].x << " " << convHull[k].y << endl;
1415 }
1416 cout << convHull[0].x << " " << convHull[0].y << endl;
1417 cout << "e" << endl;
1418 }
1419 #endif
1420
1421 #ifdef WITH_MULTISUBDIVIDE_TEST
1422 // test convex hull algorithm
1423 const double multiSubdivide_xOffset( curr_Offset );
1424 curr_Offset += 20;
1425 cout << "# multi subdivide testing" << endl
1426 << "plot [t=0:1] ";
1427 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1428 {
1429 Bezier c( someCurves[i] );
1430 Bezier c1_part1;
1431 Bezier c1_part2;
1432 Bezier c1_part3;
1433
1434 c.p0.x += multiSubdivide_xOffset;
1435 c.p1.x += multiSubdivide_xOffset;
1436 c.p2.x += multiSubdivide_xOffset;
1437 c.p3.x += multiSubdivide_xOffset;
1438
1439 const double t1( 0.1+i/(3.0*sizeof(someCurves)/sizeof(Bezier)) );
1440 const double t2( 0.9-i/(3.0*sizeof(someCurves)/sizeof(Bezier)) );
1441
1442 // subdivide at t1
1443 Impl_deCasteljauAt( c1_part1, c1_part2, c, t1 );
1444
1445 // subdivide at t2_c1. As we're working on
1446 // c1_part2 now, we have to adapt t2_c1 since
1447 // we're no longer in the original parameter
1448 // interval. This is based on the following
1449 // assumption: t2_new = (t2-t1)/(1-t1), which
1450 // relates the t2 value into the new parameter
1451 // range [0,1] of c1_part2.
1452 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
1453
1454 // subdivide at t2
1455 Impl_deCasteljauAt( c1_part3, c1_part2, c, t2 );
1456
1457 cout << " bez("
1458 << c1_part1.p0.x << ","
1459 << c1_part1.p1.x << ","
1460 << c1_part1.p2.x << ","
1461 << c1_part1.p3.x << ",t), bez("
1462 << c1_part1.p0.y+0.01 << ","
1463 << c1_part1.p1.y+0.01 << ","
1464 << c1_part1.p2.y+0.01 << ","
1465 << c1_part1.p3.y+0.01 << ",t) title \"middle " << i << "\", "
1466 << " bez("
1467 << c1_part2.p0.x << ","
1468 << c1_part2.p1.x << ","
1469 << c1_part2.p2.x << ","
1470 << c1_part2.p3.x << ",t), bez("
1471 << c1_part2.p0.y << ","
1472 << c1_part2.p1.y << ","
1473 << c1_part2.p2.y << ","
1474 << c1_part2.p3.y << ",t) title \"right " << i << "\", "
1475 << " bez("
1476 << c1_part3.p0.x << ","
1477 << c1_part3.p1.x << ","
1478 << c1_part3.p2.x << ","
1479 << c1_part3.p3.x << ",t), bez("
1480 << c1_part3.p0.y << ","
1481 << c1_part3.p1.y << ","
1482 << c1_part3.p2.y << ","
1483 << c1_part3.p3.y << ",t) title \"left " << i << "\"";
1484
1485
1486 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1487 cout << ",\\" << endl;
1488 else
1489 cout << endl;
1490 }
1491 #endif
1492
1493 #ifdef WITH_FATLINE_TEST
1494 // test fatline algorithm
1495 const double fatLine_xOffset( curr_Offset );
1496 curr_Offset += 20;
1497 cout << "# fat line testing" << endl
1498 << "plot [t=0:1] ";
1499 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1500 {
1501 Bezier c( someCurves[i] );
1502
1503 c.p0.x += fatLine_xOffset;
1504 c.p1.x += fatLine_xOffset;
1505 c.p2.x += fatLine_xOffset;
1506 c.p3.x += fatLine_xOffset;
1507
1508 FatLine line;
1509
1510 Impl_calcFatLine(line, c);
1511
1512 cout << " bez("
1513 << c.p0.x << ","
1514 << c.p1.x << ","
1515 << c.p2.x << ","
1516 << c.p3.x << ",t), bez("
1517 << c.p0.y << ","
1518 << c.p1.y << ","
1519 << c.p2.y << ","
1520 << c.p3.y << ",t) title \"bezier " << i << "\", linex("
1521 << line.a << ","
1522 << line.b << ","
1523 << line.c << ",t), liney("
1524 << line.a << ","
1525 << line.b << ","
1526 << line.c << ",t) title \"fat line (center) on " << i << "\", linex("
1527 << line.a << ","
1528 << line.b << ","
1529 << line.c-line.dMin << ",t), liney("
1530 << line.a << ","
1531 << line.b << ","
1532 << line.c-line.dMin << ",t) title \"fat line (min) on " << i << "\", linex("
1533 << line.a << ","
1534 << line.b << ","
1535 << line.c-line.dMax << ",t), liney("
1536 << line.a << ","
1537 << line.b << ","
1538 << line.c-line.dMax << ",t) title \"fat line (max) on " << i << "\"";
1539
1540 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1541 cout << ",\\" << endl;
1542 else
1543 cout << endl;
1544 }
1545 #endif
1546
1547 #ifdef WITH_CALCFOCUS_TEST
1548 // test focus curve algorithm
1549 const double focus_xOffset( curr_Offset );
1550 curr_Offset += 20;
1551 cout << "# focus line testing" << endl
1552 << "plot [t=0:1] ";
1553 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1554 {
1555 Bezier c( someCurves[i] );
1556
1557 c.p0.x += focus_xOffset;
1558 c.p1.x += focus_xOffset;
1559 c.p2.x += focus_xOffset;
1560 c.p3.x += focus_xOffset;
1561
1562 // calc focus curve
1563 Bezier focus;
1564 Impl_calcFocus(focus, c);
1565
1566 cout << " bez("
1567 << c.p0.x << ","
1568 << c.p1.x << ","
1569 << c.p2.x << ","
1570 << c.p3.x << ",t), bez("
1571 << c.p0.y << ","
1572 << c.p1.y << ","
1573 << c.p2.y << ","
1574 << c.p3.y << ",t) title \"bezier " << i << "\", bez("
1575 << focus.p0.x << ","
1576 << focus.p1.x << ","
1577 << focus.p2.x << ","
1578 << focus.p3.x << ",t), bez("
1579 << focus.p0.y << ","
1580 << focus.p1.y << ","
1581 << focus.p2.y << ","
1582 << focus.p3.y << ",t) title \"focus " << i << "\"";
1583
1584
1585 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1586 cout << ",\\" << endl;
1587 else
1588 cout << endl;
1589 }
1590 #endif
1591
1592 #ifdef WITH_SAFEPARAMBASE_TEST
1593 // test safe params base method
1594 double safeParamsBase_xOffset( curr_Offset );
1595 cout << "# safe param base method testing" << endl
1596 << "plot [t=0:1] ";
1597 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1598 {
1599 Bezier c( someCurves[i] );
1600
1601 c.p0.x += safeParamsBase_xOffset;
1602 c.p1.x += safeParamsBase_xOffset;
1603 c.p2.x += safeParamsBase_xOffset;
1604 c.p3.x += safeParamsBase_xOffset;
1605
1606 Polygon2D poly(4);
1607 poly[0] = c.p0;
1608 poly[1] = c.p1;
1609 poly[2] = c.p2;
1610 poly[3] = c.p3;
1611
1612 double t1, t2;
1613
1614 bool bRet( Impl_calcSafeParams( t1, t2, poly, 0, 1 ) );
1615
1616 Polygon2D convHull( convexHull( poly ) );
1617
1618 cout << " bez("
1619 << poly[0].x << ","
1620 << poly[1].x << ","
1621 << poly[2].x << ","
1622 << poly[3].x << ",t),bez("
1623 << poly[0].y << ","
1624 << poly[1].y << ","
1625 << poly[2].y << ","
1626 << poly[3].y << ",t), "
1627 << "t+" << safeParamsBase_xOffset << ", 0, "
1628 << "t+" << safeParamsBase_xOffset << ", 1, ";
1629 if( bRet )
1630 {
1631 cout << t1+safeParamsBase_xOffset << ", t, "
1632 << t2+safeParamsBase_xOffset << ", t, ";
1633 }
1634 cout << "'-' using ($1):($2) title \"control polygon\" with lp, "
1635 << "'-' using ($1):($2) title \"convex hull\" with lp";
1636
1637 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1638 cout << ",\\" << endl;
1639 else
1640 cout << endl;
1641
1642 safeParamsBase_xOffset += 2;
1643 }
1644
1645 safeParamsBase_xOffset = curr_Offset;
1646 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1647 {
1648 Bezier c( someCurves[i] );
1649
1650 c.p0.x += safeParamsBase_xOffset;
1651 c.p1.x += safeParamsBase_xOffset;
1652 c.p2.x += safeParamsBase_xOffset;
1653 c.p3.x += safeParamsBase_xOffset;
1654
1655 Polygon2D poly(4);
1656 poly[0] = c.p0;
1657 poly[1] = c.p1;
1658 poly[2] = c.p2;
1659 poly[3] = c.p3;
1660
1661 double t1, t2;
1662
1663 Impl_calcSafeParams( t1, t2, poly, 0, 1 );
1664
1665 Polygon2D convHull( convexHull( poly ) );
1666
1667 unsigned int k;
1668 for( k=0; k<poly.size(); ++k )
1669 {
1670 cout << poly[k].x << " " << poly[k].y << endl;
1671 }
1672 cout << poly[0].x << " " << poly[0].y << endl;
1673 cout << "e" << endl;
1674
1675 for( k=0; k<convHull.size(); ++k )
1676 {
1677 cout << convHull[k].x << " " << convHull[k].y << endl;
1678 }
1679 cout << convHull[0].x << " " << convHull[0].y << endl;
1680 cout << "e" << endl;
1681
1682 safeParamsBase_xOffset += 2;
1683 }
1684 curr_Offset += 20;
1685 #endif
1686
1687 #ifdef WITH_SAFEPARAMS_TEST
1688 // test safe parameter range algorithm
1689 const double safeParams_xOffset( curr_Offset );
1690 curr_Offset += 20;
1691 cout << "# safe param range testing" << endl
1692 << "plot [t=0.0:1.0] ";
1693 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1694 {
1695 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
1696 {
1697 Bezier c1( someCurves[i] );
1698 Bezier c2( someCurves[j] );
1699
1700 c1.p0.x += safeParams_xOffset;
1701 c1.p1.x += safeParams_xOffset;
1702 c1.p2.x += safeParams_xOffset;
1703 c1.p3.x += safeParams_xOffset;
1704 c2.p0.x += safeParams_xOffset;
1705 c2.p1.x += safeParams_xOffset;
1706 c2.p2.x += safeParams_xOffset;
1707 c2.p3.x += safeParams_xOffset;
1708
1709 double t1, t2;
1710
1711 if( Impl_calcClipRange(t1, t2, c1, c1, c2, c2) )
1712 {
1713 // clip safe ranges off c1
1714 Bezier c1_part1;
1715 Bezier c1_part2;
1716 Bezier c1_part3;
1717
1718 // subdivide at t1_c1
1719 Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 );
1720 // subdivide at t2_c1
1721 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
1722
1723 // output remaining segment (c1_part1)
1724
1725 cout << " bez("
1726 << c1.p0.x << ","
1727 << c1.p1.x << ","
1728 << c1.p2.x << ","
1729 << c1.p3.x << ",t),bez("
1730 << c1.p0.y << ","
1731 << c1.p1.y << ","
1732 << c1.p2.y << ","
1733 << c1.p3.y << ",t), bez("
1734 << c2.p0.x << ","
1735 << c2.p1.x << ","
1736 << c2.p2.x << ","
1737 << c2.p3.x << ",t),bez("
1738 << c2.p0.y << ","
1739 << c2.p1.y << ","
1740 << c2.p2.y << ","
1741 << c2.p3.y << ",t), bez("
1742 << c1_part1.p0.x << ","
1743 << c1_part1.p1.x << ","
1744 << c1_part1.p2.x << ","
1745 << c1_part1.p3.x << ",t),bez("
1746 << c1_part1.p0.y << ","
1747 << c1_part1.p1.y << ","
1748 << c1_part1.p2.y << ","
1749 << c1_part1.p3.y << ",t)";
1750
1751 if( i+2<sizeof(someCurves)/sizeof(Bezier) )
1752 cout << ",\\" << endl;
1753 else
1754 cout << endl;
1755 }
1756 }
1757 }
1758 #endif
1759
1760 #ifdef WITH_SAFEPARAM_DETAILED_TEST
1761 // test safe parameter range algorithm
1762 const double safeParams2_xOffset( curr_Offset );
1763 curr_Offset += 20;
1764 if( sizeof(someCurves)/sizeof(Bezier) > 1 )
1765 {
1766 Bezier c1( someCurves[0] );
1767 Bezier c2( someCurves[1] );
1768
1769 c1.p0.x += safeParams2_xOffset;
1770 c1.p1.x += safeParams2_xOffset;
1771 c1.p2.x += safeParams2_xOffset;
1772 c1.p3.x += safeParams2_xOffset;
1773 c2.p0.x += safeParams2_xOffset;
1774 c2.p1.x += safeParams2_xOffset;
1775 c2.p2.x += safeParams2_xOffset;
1776 c2.p3.x += safeParams2_xOffset;
1777
1778 double t1, t2;
1779
1780 // output happens here
1781 Impl_calcClipRange(t1, t2, c1, c1, c2, c2);
1782 }
1783 #endif
1784
1785 #ifdef WITH_SAFEFOCUSPARAM_TEST
1786 // test safe parameter range from focus algorithm
1787 const double safeParamsFocus_xOffset( curr_Offset );
1788 curr_Offset += 20;
1789 cout << "# safe param range from focus testing" << endl
1790 << "plot [t=0.0:1.0] ";
1791 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1792 {
1793 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
1794 {
1795 Bezier c1( someCurves[i] );
1796 Bezier c2( someCurves[j] );
1797
1798 c1.p0.x += safeParamsFocus_xOffset;
1799 c1.p1.x += safeParamsFocus_xOffset;
1800 c1.p2.x += safeParamsFocus_xOffset;
1801 c1.p3.x += safeParamsFocus_xOffset;
1802 c2.p0.x += safeParamsFocus_xOffset;
1803 c2.p1.x += safeParamsFocus_xOffset;
1804 c2.p2.x += safeParamsFocus_xOffset;
1805 c2.p3.x += safeParamsFocus_xOffset;
1806
1807 double t1, t2;
1808
1809 Bezier focus;
1810 #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
1811 #if 0
1812 {
1813 // clip safe ranges off c1_orig
1814 Bezier c1_part1;
1815 Bezier c1_part2;
1816 Bezier c1_part3;
1817
1818 // subdivide at t1_c1
1819 Impl_deCasteljauAt( c1_part1, c1_part2, c2, 0.30204 );
1820
1821 // subdivide at t2_c1. As we're working on
1822 // c1_part2 now, we have to adapt t2_c1 since
1823 // we're no longer in the original parameter
1824 // interval. This is based on the following
1825 // assumption: t2_new = (t2-t1)/(1-t1), which
1826 // relates the t2 value into the new parameter
1827 // range [0,1] of c1_part2.
1828 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (0.57151-0.30204)/(1.0-0.30204) );
1829
1830 c2 = c1_part1;
1831 Impl_calcFocus( focus, c2 );
1832 }
1833 #else
1834 Impl_calcFocus( focus, c2 );
1835 #endif
1836 #else
1837 focus = c2;
1838 #endif
1839 // determine safe range on c1
1840 bool bRet( Impl_calcSafeParams_focus( t1, t2,
1841 c1, focus ) );
1842
1843 cerr << "t1: " << t1 << ", t2: " << t2 << endl;
1844
1845 // clip safe ranges off c1
1846 Bezier c1_part1;
1847 Bezier c1_part2;
1848 Bezier c1_part3;
1849
1850 // subdivide at t1_c1
1851 Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 );
1852 // subdivide at t2_c1
1853 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
1854
1855 // output remaining segment (c1_part1)
1856
1857 cout << " bez("
1858 << c1.p0.x << ","
1859 << c1.p1.x << ","
1860 << c1.p2.x << ","
1861 << c1.p3.x << ",t),bez("
1862 << c1.p0.y << ","
1863 << c1.p1.y << ","
1864 << c1.p2.y << ","
1865 << c1.p3.y << ",t) title \"c1\", "
1866 #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
1867 << "bez("
1868 << c2.p0.x << ","
1869 << c2.p1.x << ","
1870 << c2.p2.x << ","
1871 << c2.p3.x << ",t),bez("
1872 << c2.p0.y << ","
1873 << c2.p1.y << ","
1874 << c2.p2.y << ","
1875 << c2.p3.y << ",t) title \"c2\", "
1876 << "bez("
1877 << focus.p0.x << ","
1878 << focus.p1.x << ","
1879 << focus.p2.x << ","
1880 << focus.p3.x << ",t),bez("
1881 << focus.p0.y << ","
1882 << focus.p1.y << ","
1883 << focus.p2.y << ","
1884 << focus.p3.y << ",t) title \"focus\"";
1885 #else
1886 << "bez("
1887 << c2.p0.x << ","
1888 << c2.p1.x << ","
1889 << c2.p2.x << ","
1890 << c2.p3.x << ",t),bez("
1891 << c2.p0.y << ","
1892 << c2.p1.y << ","
1893 << c2.p2.y << ","
1894 << c2.p3.y << ",t) title \"focus\"";
1895 #endif
1896 if( bRet )
1897 {
1898 cout << ", bez("
1899 << c1_part1.p0.x << ","
1900 << c1_part1.p1.x << ","
1901 << c1_part1.p2.x << ","
1902 << c1_part1.p3.x << ",t),bez("
1903 << c1_part1.p0.y+0.01 << ","
1904 << c1_part1.p1.y+0.01 << ","
1905 << c1_part1.p2.y+0.01 << ","
1906 << c1_part1.p3.y+0.01 << ",t) title \"part\"";
1907 }
1908
1909 if( i+2<sizeof(someCurves)/sizeof(Bezier) )
1910 cout << ",\\" << endl;
1911 else
1912 cout << endl;
1913 }
1914 }
1915 #endif
1916
1917 #ifdef WITH_SAFEFOCUSPARAM_DETAILED_TEST
1918 // test safe parameter range algorithm
1919 const double safeParams3_xOffset( curr_Offset );
1920 curr_Offset += 20;
1921 if( sizeof(someCurves)/sizeof(Bezier) > 1 )
1922 {
1923 Bezier c1( someCurves[0] );
1924 Bezier c2( someCurves[1] );
1925
1926 c1.p0.x += safeParams3_xOffset;
1927 c1.p1.x += safeParams3_xOffset;
1928 c1.p2.x += safeParams3_xOffset;
1929 c1.p3.x += safeParams3_xOffset;
1930 c2.p0.x += safeParams3_xOffset;
1931 c2.p1.x += safeParams3_xOffset;
1932 c2.p2.x += safeParams3_xOffset;
1933 c2.p3.x += safeParams3_xOffset;
1934
1935 double t1, t2;
1936
1937 Bezier focus;
1938 #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
1939 Impl_calcFocus( focus, c2 );
1940 #else
1941 focus = c2;
1942 #endif
1943
1944 // determine safe range on c1, output happens here
1945 Impl_calcSafeParams_focus( t1, t2,
1946 c1, focus );
1947 }
1948 #endif
1949
1950 #ifdef WITH_BEZIERCLIP_TEST
1951 ::std::vector< ::std::pair<double, double> > result;
1952 ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > > ii(result);
1953
1954 // test full bezier clipping
1955 const double bezierClip_xOffset( curr_Offset );
1956 curr_Offset += 20;
1957 cout << endl << endl << "# bezier clip testing" << endl
1958 << "plot [t=0:1] ";
1959 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1960 {
1961 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
1962 {
1963 Bezier c1( someCurves[i] );
1964 Bezier c2( someCurves[j] );
1965
1966 c1.p0.x += bezierClip_xOffset;
1967 c1.p1.x += bezierClip_xOffset;
1968 c1.p2.x += bezierClip_xOffset;
1969 c1.p3.x += bezierClip_xOffset;
1970 c2.p0.x += bezierClip_xOffset;
1971 c2.p1.x += bezierClip_xOffset;
1972 c2.p2.x += bezierClip_xOffset;
1973 c2.p3.x += bezierClip_xOffset;
1974
1975 cout << " bez("
1976 << c1.p0.x << ","
1977 << c1.p1.x << ","
1978 << c1.p2.x << ","
1979 << c1.p3.x << ",t),bez("
1980 << c1.p0.y << ","
1981 << c1.p1.y << ","
1982 << c1.p2.y << ","
1983 << c1.p3.y << ",t), bez("
1984 << c2.p0.x << ","
1985 << c2.p1.x << ","
1986 << c2.p2.x << ","
1987 << c2.p3.x << ",t),bez("
1988 << c2.p0.y << ","
1989 << c2.p1.y << ","
1990 << c2.p2.y << ","
1991 << c2.p3.y << ",t), '-' using (bez("
1992 << c1.p0.x << ","
1993 << c1.p1.x << ","
1994 << c1.p2.x << ","
1995 << c1.p3.x
1996 << ",$1)):(bez("
1997 << c1.p0.y << ","
1998 << c1.p1.y << ","
1999 << c1.p2.y << ","
2000 << c1.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << i << ")\", "
2001 << " '-' using (bez("
2002 << c2.p0.x << ","
2003 << c2.p1.x << ","
2004 << c2.p2.x << ","
2005 << c2.p3.x
2006 << ",$1)):(bez("
2007 << c2.p0.y << ","
2008 << c2.p1.y << ","
2009 << c2.p2.y << ","
2010 << c2.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << j << ")\"";
2011
2012 if( i+2<sizeof(someCurves)/sizeof(Bezier) )
2013 cout << ",\\" << endl;
2014 else
2015 cout << endl;
2016 }
2017 }
2018 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
2019 {
2020 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
2021 {
2022 result.clear();
2023 Bezier c1( someCurves[i] );
2024 Bezier c2( someCurves[j] );
2025
2026 c1.p0.x += bezierClip_xOffset;
2027 c1.p1.x += bezierClip_xOffset;
2028 c1.p2.x += bezierClip_xOffset;
2029 c1.p3.x += bezierClip_xOffset;
2030 c2.p0.x += bezierClip_xOffset;
2031 c2.p1.x += bezierClip_xOffset;
2032 c2.p2.x += bezierClip_xOffset;
2033 c2.p3.x += bezierClip_xOffset;
2034
2035 clipBezier( ii, 0.00001, c1, c2 );
2036
2037 for( k=0; k<result.size(); ++k )
2038 {
2039 cout << result[k].first << endl;
2040 }
2041 cout << "e" << endl;
2042
2043 for( k=0; k<result.size(); ++k )
2044 {
2045 cout << result[k].second << endl;
2046 }
2047 cout << "e" << endl;
2048 }
2049 }
2050 #endif
2051
2052 return 0;
2053 }
2054