update dev300-m58
[ooovba.git] / basegfx / source / workbench / bezierclip.cxx
blobc5e49f6651edfd659788a380ba295406ba764733
1 /*************************************************************************
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * Copyright 2008 by Sun Microsystems, Inc.
7 * OpenOffice.org - a multi-platform office productivity suite
9 * $RCSfile: bezierclip.cxx,v $
10 * $Revision: 1.6 $
12 * This file is part of OpenOffice.org.
14 * OpenOffice.org is free software: you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License version 3
16 * only, as published by the Free Software Foundation.
18 * OpenOffice.org is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU Lesser General Public License version 3 for more details
22 * (a copy is included in the LICENSE file that accompanied this code).
24 * You should have received a copy of the GNU Lesser General Public License
25 * version 3 along with OpenOffice.org. If not, see
26 * <http://www.openoffice.org/license.html>
27 * for a copy of the LGPLv3 License.
29 ************************************************************************/
31 // MARKER(update_precomp.py): autogen include statement, do not remove
32 #include "precompiled_basegfx.hxx"
34 #include <algorithm>
35 #include <iterator>
36 #include <vector>
37 #include <utility>
39 #include <math.h>
41 #include "bezierclip.hxx"
42 #include "gauss.hxx"
46 // what to test
47 #define WITH_ASSERTIONS
48 //#define WITH_CONVEXHULL_TEST
49 //#define WITH_MULTISUBDIVIDE_TEST
50 //#define WITH_FATLINE_TEST
51 //#define WITH_CALCFOCUS_TEST
52 //#define WITH_SAFEPARAMBASE_TEST
53 //#define WITH_SAFEPARAMS_TEST
54 //#define WITH_SAFEPARAM_DETAILED_TEST
55 //#define WITH_SAFEFOCUSPARAM_CALCFOCUS
56 //#define WITH_SAFEFOCUSPARAM_TEST
57 //#define WITH_SAFEFOCUSPARAM_DETAILED_TEST
58 #define WITH_BEZIERCLIP_TEST
62 // -----------------------------------------------------------------------------
64 /* Implementation of the so-called 'Fat-Line Bezier Clipping Algorithm' by Sederberg et al.
66 * Actual reference is: T. W. Sederberg and T Nishita: Curve
67 * intersection using Bezier clipping. In Computer Aided Design, 22
68 * (9), 1990, pp. 538--549
71 // -----------------------------------------------------------------------------
73 /* Misc helper
74 * ===========
76 int fallFac( int n, int k )
78 #ifdef WITH_ASSERTIONS
79 assert(n>=k); // "For factorials, n must be greater or equal k"
80 assert(n>=0); // "For factorials, n must be positive"
81 assert(k>=0); // "For factorials, k must be positive"
82 #endif
84 int res( 1 );
86 while( k-- && n ) res *= n--;
88 return res;
91 // -----------------------------------------------------------------------------
93 int fac( int n )
95 return fallFac(n, n);
98 // -----------------------------------------------------------------------------
100 /* Bezier fat line clipping part
101 * =============================
104 // -----------------------------------------------------------------------------
106 void Impl_calcFatLine( FatLine& line, const Bezier& c )
108 // Prepare normalized implicit line
109 // ================================
111 // calculate vector orthogonal to p1-p4:
112 line.a = -(c.p0.y - c.p3.y);
113 line.b = (c.p0.x - c.p3.x);
115 // normalize
116 const double len( sqrt( line.a*line.a + line.b*line.b ) );
117 if( !tolZero(len) )
119 line.a /= len;
120 line.b /= len;
123 line.c = -(line.a*c.p0.x + line.b*c.p0.y);
126 // Determine bounding fat line from it
127 // ===================================
129 // calc control point distances
130 const double dP2( calcLineDistance(line.a, line.b, line.c, c.p1.x, c.p1.y ) );
131 const double dP3( calcLineDistance(line.a, line.b, line.c, c.p2.x, c.p2.y ) );
133 // calc approximate bounding lines to curve (tight bounds are
134 // possible here, but more expensive to calculate and thus not
135 // worth the overhead)
136 if( dP2 * dP3 > 0.0 )
138 line.dMin = 3.0/4.0 * ::std::min(0.0, ::std::min(dP2, dP3));
139 line.dMax = 3.0/4.0 * ::std::max(0.0, ::std::max(dP2, dP3));
141 else
143 line.dMin = 4.0/9.0 * ::std::min(0.0, ::std::min(dP2, dP3));
144 line.dMax = 4.0/9.0 * ::std::max(0.0, ::std::max(dP2, dP3));
148 void Impl_calcBounds( Point2D& leftTop,
149 Point2D& rightBottom,
150 const Bezier& c1 )
152 leftTop.x = ::std::min( c1.p0.x, ::std::min( c1.p1.x, ::std::min( c1.p2.x, c1.p3.x ) ) );
153 leftTop.y = ::std::min( c1.p0.y, ::std::min( c1.p1.y, ::std::min( c1.p2.y, c1.p3.y ) ) );
154 rightBottom.x = ::std::max( c1.p0.x, ::std::max( c1.p1.x, ::std::max( c1.p2.x, c1.p3.x ) ) );
155 rightBottom.y = ::std::max( c1.p0.y, ::std::max( c1.p1.y, ::std::max( c1.p2.y, c1.p3.y ) ) );
158 bool Impl_doBBoxIntersect( const Bezier& c1,
159 const Bezier& c2 )
161 // calc rectangular boxes from c1 and c2
162 Point2D lt1;
163 Point2D rb1;
164 Point2D lt2;
165 Point2D rb2;
167 Impl_calcBounds( lt1, rb1, c1 );
168 Impl_calcBounds( lt2, rb2, c2 );
170 if( ::std::min(rb1.x, rb2.x) < ::std::max(lt1.x, lt2.x) ||
171 ::std::min(rb1.y, rb2.y) < ::std::max(lt1.y, lt2.y) )
173 return false;
175 else
177 return true;
181 /* calculates two t's for the given bernstein control polygon: the first is
182 * the intersection of the min value line with the convex hull from
183 * the left, the second is the intersection of the max value line with
184 * the convex hull from the right.
186 bool Impl_calcSafeParams( double& t1,
187 double& t2,
188 const Polygon2D& rPoly,
189 double lowerYBound,
190 double upperYBound )
192 // need the convex hull of the control polygon, as this is
193 // guaranteed to completely bound the curve
194 Polygon2D convHull( convexHull(rPoly) );
196 // init min and max buffers
197 t1 = 0.0 ;
198 double currLowerT( 1.0 );
200 t2 = 1.0;
201 double currHigherT( 0.0 );
203 if( convHull.size() <= 1 )
204 return false; // only one point? Then we're done with clipping
206 /* now, clip against lower and higher bounds */
207 Point2D p0;
208 Point2D p1;
210 bool bIntersection( false );
212 for( Polygon2D::size_type i=0; i<convHull.size(); ++i )
214 // have to check against convHull.size() segments, as the
215 // convex hull is, by definition, closed. Thus, for the
216 // last point, we take the first point as partner.
217 if( i+1 == convHull.size() )
219 // close the polygon
220 p0 = convHull[i];
221 p1 = convHull[0];
223 else
225 p0 = convHull[i];
226 p1 = convHull[i+1];
229 // is the segment in question within or crossing the
230 // horizontal band spanned by lowerYBound and upperYBound? If
231 // not, we've got no intersection. If yes, we maybe don't have
232 // an intersection, but we've got to update the permissible
233 // range, nevertheless. This is because inside lying segments
234 // leads to full range forbidden.
235 if( (tolLessEqual(p0.y, upperYBound) || tolLessEqual(p1.y, upperYBound)) &&
236 (tolGreaterEqual(p0.y, lowerYBound) || tolGreaterEqual(p1.y, lowerYBound)) )
238 // calc intersection of convex hull segment with
239 // one of the horizontal bounds lines
240 const double r_x( p1.x - p0.x );
241 const double r_y( p1.y - p0.y );
243 if( tolZero(r_y) )
245 // r_y is virtually zero, thus we've got a horizontal
246 // line. Now check whether we maybe coincide with lower or
247 // upper horizonal bound line.
248 if( tolEqual(p0.y, lowerYBound) ||
249 tolEqual(p0.y, upperYBound) )
251 // yes, simulate intersection then
252 currLowerT = ::std::min(currLowerT, ::std::min(p0.x, p1.x));
253 currHigherT = ::std::max(currHigherT, ::std::max(p0.x, p1.x));
256 else
258 // check against lower and higher bounds
259 // =====================================
261 // calc intersection with horizontal dMin line
262 const double currTLow( (lowerYBound - p0.y) * r_x / r_y + p0.x );
264 // calc intersection with horizontal dMax line
265 const double currTHigh( (upperYBound - p0.y) * r_x / r_y + p0.x );
267 currLowerT = ::std::min(currLowerT, ::std::min(currTLow, currTHigh));
268 currHigherT = ::std::max(currHigherT, ::std::max(currTLow, currTHigh));
271 // set flag that at least one segment is contained or
272 // intersects given horizontal band.
273 bIntersection = true;
277 #ifndef WITH_SAFEPARAMBASE_TEST
278 // limit intersections found to permissible t parameter range
279 t1 = ::std::max(0.0, currLowerT);
280 t2 = ::std::min(1.0, currHigherT);
281 #endif
283 return bIntersection;
287 /* calculates two t's for the given bernstein polynomial: the first is
288 * the intersection of the min value line with the convex hull from
289 * the left, the second is the intersection of the max value line with
290 * the convex hull from the right.
292 * The polynomial coefficients c0 to c3 given to this method
293 * must correspond to t values of 0, 1/3, 2/3 and 1, respectively.
295 bool Impl_calcSafeParams_clip( double& t1,
296 double& t2,
297 const FatLine& bounds,
298 double c0,
299 double c1,
300 double c2,
301 double c3 )
303 /* first of all, determine convex hull of c0-c3 */
304 Polygon2D poly(4);
305 poly[0] = Point2D(0, c0);
306 poly[1] = Point2D(1.0/3.0, c1);
307 poly[2] = Point2D(2.0/3.0, c2);
308 poly[3] = Point2D(1, c3);
310 #ifndef WITH_SAFEPARAM_DETAILED_TEST
312 return Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax );
314 #else
315 bool bRet( Impl_calcSafeParams( t1, t2, poly, bounds.dMin, bounds.dMax ) );
317 Polygon2D convHull( convexHull( poly ) );
319 cout << "# convex hull testing" << endl
320 << "plot [t=0:1] ";
321 cout << " bez("
322 << poly[0].x << ","
323 << poly[1].x << ","
324 << poly[2].x << ","
325 << poly[3].x << ",t),bez("
326 << poly[0].y << ","
327 << poly[1].y << ","
328 << poly[2].y << ","
329 << poly[3].y << ",t), "
330 << "t, " << bounds.dMin << ", "
331 << "t, " << bounds.dMax << ", "
332 << t1 << ", t, "
333 << t2 << ", t, "
334 << "'-' using ($1):($2) title \"control polygon\" with lp, "
335 << "'-' using ($1):($2) title \"convex hull\" with lp" << endl;
337 unsigned int k;
338 for( k=0; k<poly.size(); ++k )
340 cout << poly[k].x << " " << poly[k].y << endl;
342 cout << poly[0].x << " " << poly[0].y << endl;
343 cout << "e" << endl;
345 for( k=0; k<convHull.size(); ++k )
347 cout << convHull[k].x << " " << convHull[k].y << endl;
349 cout << convHull[0].x << " " << convHull[0].y << endl;
350 cout << "e" << endl;
352 return bRet;
353 #endif
356 // -----------------------------------------------------------------------------
358 void Impl_deCasteljauAt( Bezier& part1,
359 Bezier& part2,
360 const Bezier& input,
361 double t )
363 // deCasteljau bezier arc, scheme is:
365 // First row is C_0^n,C_1^n,...,C_n^n
366 // Second row is P_1^n,...,P_n^n
367 // etc.
368 // with P_k^r = (1 - x_s)P_{k-1}^{r-1} + x_s P_k{r-1}
370 // this results in:
372 // P1 P2 P3 P4
373 // L1 P2 P3 R4
374 // L2 H R3
375 // L3 R2
376 // L4/R1
377 if( tolZero(t) )
379 // t is zero -> part2 is input curve, part1 is empty (input.p0, that is)
380 part1.p0.x = part1.p1.x = part1.p2.x = part1.p3.x = input.p0.x;
381 part1.p0.y = part1.p1.y = part1.p2.y = part1.p3.y = input.p0.y;
382 part2 = input;
384 else if( tolEqual(t, 1.0) )
386 // t is one -> part1 is input curve, part2 is empty (input.p3, that is)
387 part1 = input;
388 part2.p0.x = part2.p1.x = part2.p2.x = part2.p3.x = input.p3.x;
389 part2.p0.y = part2.p1.y = part2.p2.y = part2.p3.y = input.p3.y;
391 else
393 part1.p0.x = input.p0.x; part1.p0.y = input.p0.y;
394 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;
395 const double Hx ( (1.0 - t)*input.p1.x + t*input.p2.x ), Hy ( (1.0 - t)*input.p1.y + t*input.p2.y );
396 part1.p2.x = (1.0 - t)*part1.p1.x + t*Hx; part1.p2.y = (1.0 - t)*part1.p1.y + t*Hy;
397 part2.p3.x = input.p3.x; part2.p3.y = input.p3.y;
398 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;
399 part2.p1.x = (1.0 - t)*Hx + t*part2.p2.x; part2.p1.y = (1.0 - t)*Hy + t*part2.p2.y;
400 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;
401 part1.p3.x = part2.p0.x; part1.p3.y = part2.p0.y;
405 // -----------------------------------------------------------------------------
407 void printCurvesWithSafeRange( const Bezier& c1, const Bezier& c2, double t1_c1, double t2_c1,
408 const Bezier& c2_part, const FatLine& bounds_c2 )
410 static int offset = 0;
412 cout << "# safe param range testing" << endl
413 << "plot [t=0.0:1.0] ";
415 // clip safe ranges off c1
416 Bezier c1_part1;
417 Bezier c1_part2;
418 Bezier c1_part3;
420 // subdivide at t1_c1
421 Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1_c1 );
422 // subdivide at t2_c1
423 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, t2_c1 );
425 // output remaining segment (c1_part1)
427 cout << "bez("
428 << c1.p0.x+offset << ","
429 << c1.p1.x+offset << ","
430 << c1.p2.x+offset << ","
431 << c1.p3.x+offset << ",t),bez("
432 << c1.p0.y << ","
433 << c1.p1.y << ","
434 << c1.p2.y << ","
435 << c1.p3.y << ",t), bez("
436 << c2.p0.x+offset << ","
437 << c2.p1.x+offset << ","
438 << c2.p2.x+offset << ","
439 << c2.p3.x+offset << ",t),bez("
440 << c2.p0.y << ","
441 << c2.p1.y << ","
442 << c2.p2.y << ","
443 << c2.p3.y << ",t), "
444 #if 1
445 << "bez("
446 << c1_part1.p0.x+offset << ","
447 << c1_part1.p1.x+offset << ","
448 << c1_part1.p2.x+offset << ","
449 << c1_part1.p3.x+offset << ",t),bez("
450 << c1_part1.p0.y << ","
451 << c1_part1.p1.y << ","
452 << c1_part1.p2.y << ","
453 << c1_part1.p3.y << ",t), "
454 #endif
455 #if 1
456 << "bez("
457 << c2_part.p0.x+offset << ","
458 << c2_part.p1.x+offset << ","
459 << c2_part.p2.x+offset << ","
460 << c2_part.p3.x+offset << ",t),bez("
461 << c2_part.p0.y << ","
462 << c2_part.p1.y << ","
463 << c2_part.p2.y << ","
464 << c2_part.p3.y << ",t), "
465 #endif
466 << "linex("
467 << bounds_c2.a << ","
468 << bounds_c2.b << ","
469 << bounds_c2.c << ",t)+" << offset << ", liney("
470 << bounds_c2.a << ","
471 << bounds_c2.b << ","
472 << bounds_c2.c << ",t) title \"fat line (center)\", linex("
473 << bounds_c2.a << ","
474 << bounds_c2.b << ","
475 << bounds_c2.c-bounds_c2.dMin << ",t)+" << offset << ", liney("
476 << bounds_c2.a << ","
477 << bounds_c2.b << ","
478 << bounds_c2.c-bounds_c2.dMin << ",t) title \"fat line (min) \", linex("
479 << bounds_c2.a << ","
480 << bounds_c2.b << ","
481 << bounds_c2.c-bounds_c2.dMax << ",t)+" << offset << ", liney("
482 << bounds_c2.a << ","
483 << bounds_c2.b << ","
484 << bounds_c2.c-bounds_c2.dMax << ",t) title \"fat line (max) \"" << endl;
486 offset += 1;
489 // -----------------------------------------------------------------------------
491 void printResultWithFinalCurves( const Bezier& c1, const Bezier& c1_part,
492 const Bezier& c2, const Bezier& c2_part,
493 double t1_c1, double t2_c1 )
495 static int offset = 0;
497 cout << "# final result" << endl
498 << "plot [t=0.0:1.0] ";
500 cout << "bez("
501 << c1.p0.x+offset << ","
502 << c1.p1.x+offset << ","
503 << c1.p2.x+offset << ","
504 << c1.p3.x+offset << ",t),bez("
505 << c1.p0.y << ","
506 << c1.p1.y << ","
507 << c1.p2.y << ","
508 << c1.p3.y << ",t), bez("
509 << c1_part.p0.x+offset << ","
510 << c1_part.p1.x+offset << ","
511 << c1_part.p2.x+offset << ","
512 << c1_part.p3.x+offset << ",t),bez("
513 << c1_part.p0.y << ","
514 << c1_part.p1.y << ","
515 << c1_part.p2.y << ","
516 << c1_part.p3.y << ",t), "
517 << " pointmarkx(bez("
518 << c1.p0.x+offset << ","
519 << c1.p1.x+offset << ","
520 << c1.p2.x+offset << ","
521 << c1.p3.x+offset << ","
522 << t1_c1 << "),t), "
523 << " pointmarky(bez("
524 << c1.p0.y << ","
525 << c1.p1.y << ","
526 << c1.p2.y << ","
527 << c1.p3.y << ","
528 << t1_c1 << "),t), "
529 << " pointmarkx(bez("
530 << c1.p0.x+offset << ","
531 << c1.p1.x+offset << ","
532 << c1.p2.x+offset << ","
533 << c1.p3.x+offset << ","
534 << t2_c1 << "),t), "
535 << " pointmarky(bez("
536 << c1.p0.y << ","
537 << c1.p1.y << ","
538 << c1.p2.y << ","
539 << c1.p3.y << ","
540 << t2_c1 << "),t), "
542 << "bez("
543 << c2.p0.x+offset << ","
544 << c2.p1.x+offset << ","
545 << c2.p2.x+offset << ","
546 << c2.p3.x+offset << ",t),bez("
547 << c2.p0.y << ","
548 << c2.p1.y << ","
549 << c2.p2.y << ","
550 << c2.p3.y << ",t), "
551 << "bez("
552 << c2_part.p0.x+offset << ","
553 << c2_part.p1.x+offset << ","
554 << c2_part.p2.x+offset << ","
555 << c2_part.p3.x+offset << ",t),bez("
556 << c2_part.p0.y << ","
557 << c2_part.p1.y << ","
558 << c2_part.p2.y << ","
559 << c2_part.p3.y << ",t)" << endl;
561 offset += 1;
564 // -----------------------------------------------------------------------------
566 /** determine parameter ranges [0,t1) and (t2,1] on c1, where c1 is guaranteed to lie outside c2.
567 Returns false, if the two curves don't even intersect.
569 @param t1
570 Range [0,t1) on c1 is guaranteed to lie outside c2
572 @param t2
573 Range (t2,1] on c1 is guaranteed to lie outside c2
575 @param c1_orig
576 Original curve c1
578 @param c1_part
579 Subdivided current part of c1
581 @param c2_orig
582 Original curve c2
584 @param c2_part
585 Subdivided current part of c2
587 bool Impl_calcClipRange( double& t1,
588 double& t2,
589 const Bezier& c1_orig,
590 const Bezier& c1_part,
591 const Bezier& c2_orig,
592 const Bezier& c2_part )
594 // TODO: Maybe also check fat line orthogonal to P0P3, having P0
595 // and P3 as the extremal points
597 if( Impl_doBBoxIntersect(c1_part, c2_part) )
599 // Calculate fat lines around c1
600 FatLine bounds_c2;
602 // must use the subdivided version of c2, since the fat line
603 // algorithm works implicitely with the convex hull bounding
604 // box.
605 Impl_calcFatLine(bounds_c2, c2_part);
607 // determine clip positions on c2. Can use original c1 (which
608 // is necessary anyway, to get the t's on the original curve),
609 // since the distance calculations work directly in the
610 // Bernstein polynom parameter domain.
611 if( Impl_calcSafeParams_clip( t1, t2, bounds_c2,
612 calcLineDistance( bounds_c2.a,
613 bounds_c2.b,
614 bounds_c2.c,
615 c1_orig.p0.x,
616 c1_orig.p0.y ),
617 calcLineDistance( bounds_c2.a,
618 bounds_c2.b,
619 bounds_c2.c,
620 c1_orig.p1.x,
621 c1_orig.p1.y ),
622 calcLineDistance( bounds_c2.a,
623 bounds_c2.b,
624 bounds_c2.c,
625 c1_orig.p2.x,
626 c1_orig.p2.y ),
627 calcLineDistance( bounds_c2.a,
628 bounds_c2.b,
629 bounds_c2.c,
630 c1_orig.p3.x,
631 c1_orig.p3.y ) ) )
633 //printCurvesWithSafeRange(c1_orig, c2_orig, t1, t2, c2_part, bounds_c2);
635 // they do intersect
636 return true;
640 // they don't intersect: nothing to do
641 return false;
644 // -----------------------------------------------------------------------------
646 /* Tangent intersection part
647 * =========================
650 // -----------------------------------------------------------------------------
652 void Impl_calcFocus( Bezier& res, const Bezier& c )
654 // arbitrary small value, for now
655 // TODO: find meaningful value
656 const double minPivotValue( 1.0e-20 );
658 Point2D::value_type fMatrix[6];
659 Point2D::value_type fRes[2];
661 // calc new curve from hodograph, c and linear blend
663 // Coefficients for derivative of c are (C_i=n(C_{i+1} - C_i)):
665 // 3(P1 - P0), 3(P2 - P1), 3(P3 - P2) (bezier curve of degree 2)
667 // The hodograph is then (bezier curve of 2nd degree is P0(1-t)^2 + 2P1(1-t)t + P2t^2):
669 // 3(P1 - P0)(1-t)^2 + 6(P2 - P1)(1-t)t + 3(P3 - P2)t^2
671 // rotate by 90 degrees: x=-y, y=x and you get the normal vector function N(t):
673 // 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)
674 // 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
676 // Now, the focus curve is defined to be F(t)=P(t) + c(t)N(t),
677 // where P(t) is the original curve, and c(t)=c0(1-t) + c1 t
679 // This results in the following expression for F(t):
681 // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
682 // (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)
684 // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1.t)t^2 + P3.y t^3 +
685 // (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)
687 // As a heuristic, we set F(0)=F(1) (thus, the curve is closed and _tends_ to be small):
689 // For F(0), the following results:
691 // x(0) = P0.x - c0 3(P1.y - P0.y)
692 // y(0) = P0.y + c0 3(P1.x - P0.x)
694 // For F(1), the following results:
696 // x(1) = P3.x - c1 3(P3.y - P2.y)
697 // y(1) = P3.y + c1 3(P3.x - P2.x)
699 // Reorder, collect and substitute into F(0)=F(1):
701 // P0.x - c0 3(P1.y - P0.y) = P3.x - c1 3(P3.y - P2.y)
702 // P0.y + c0 3(P1.x - P0.x) = P3.y + c1 3(P3.x - P2.x)
704 // which yields
706 // (P0.y - P1.y)c0 + (P3.y - P2.y)c1 = (P3.x - P0.x)/3
707 // (P1.x - P0.x)c0 + (P2.x - P3.x)c1 = (P3.y - P0.y)/3
710 // so, this is what we calculate here (determine c0 and c1):
711 fMatrix[0] = c.p1.x - c.p0.x;
712 fMatrix[1] = c.p2.x - c.p3.x;
713 fMatrix[2] = (c.p3.y - c.p0.y)/3.0;
714 fMatrix[3] = c.p0.y - c.p1.y;
715 fMatrix[4] = c.p3.y - c.p2.y;
716 fMatrix[5] = (c.p3.x - c.p0.x)/3.0;
718 // TODO: determine meaningful value for
719 if( !solve(fMatrix, 2, 3, fRes, minPivotValue) )
721 // TODO: generate meaningful values here
722 // singular or nearly singular system -- use arbitrary
723 // values for res
724 fRes[0] = 0.0;
725 fRes[1] = 1.0;
727 cerr << "Matrix singular!" << endl;
730 // now, the reordered and per-coefficient collected focus curve is
731 // the following third degree bezier curve F(t):
733 // x(t) = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
734 // (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)
735 // = P0.x (1-t)^3 + 3 P1.x (1-t)^2t + 3 P2.x (1.t)t^2 + P3.x t^3 -
736 // (3c0P1.y(1-t)^3 - 3c0P0.y(1-t)^3 + 6c0P2.y(1-t)^2t - 6c0P1.y(1-t)^2t +
737 // 3c0P3.y(1-t)t^2 - 3c0P2.y(1-t)t^2 +
738 // 3c1P1.y(1-t)^2t - 3c1P0.y(1-t)^2t + 6c1P2.y(1-t)t^2 - 6c1P1.y(1-t)t^2 +
739 // 3c1P3.yt^3 - 3c1P2.yt^3)
740 // = (P0.x - 3 c0 P1.y + 3 c0 P0.y)(1-t)^3 +
741 // 3(P1.x - c1 P1.y + c1 P0.y - 2 c0 P2.y + 2 c0 P1.y)(1-t)^2t +
742 // 3(P2.x - 2 c1 P2.y + 2 c1 P1.y - c0 P3.y + c0 P2.y)(1-t)t^2 +
743 // (P3.x - 3 c1 P3.y + 3 c1 P2.y)t^3
744 // = (P0.x - 3 c0(P1.y - P0.y))(1-t)^3 +
745 // 3(P1.x - c1(P1.y - P0.y) - 2c0(P2.y - P1.y))(1-t)^2t +
746 // 3(P2.x - 2 c1(P2.y - P1.y) - c0(P3.y - P2.y))(1-t)t^2 +
747 // (P3.x - 3 c1(P3.y - P2.y))t^3
749 // y(t) = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 +
750 // (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)
751 // = P0.y (1-t)^3 + 3 P1.y (1-t)^2t + 3 P2.y (1-t)t^2 + P3.y t^3 +
752 // 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 +
753 // 3c1(P1.x - P0.x)(1-t)^2t + 6c1(P2.x - P1.x)(1-t)t^2 + 3c1(P3.x - P2.x)t^3
754 // = (P0.y + 3 c0 (P1.x - P0.x))(1-t)^3 +
755 // 3(P1.y + 2 c0 (P2.x - P1.x) + c1 (P1.x - P0.x))(1-t)^2t +
756 // 3(P2.y + c0 (P3.x - P2.x) + 2 c1 (P2.x - P1.x))(1-t)t^2 +
757 // (P3.y + 3 c1 (P3.x - P2.x))t^3
759 // Therefore, the coefficients F0 to F3 of the focus curve are:
761 // F0.x = (P0.x - 3 c0(P1.y - P0.y)) F0.y = (P0.y + 3 c0 (P1.x - P0.x))
762 // 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))
763 // 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))
764 // F3.x = (P3.x - 3 c1(P3.y - P2.y)) F3.y = (P3.y + 3 c1 (P3.x - P2.x))
766 res.p0.x = c.p0.x - 3*fRes[0]*(c.p1.y - c.p0.y);
767 res.p1.x = c.p1.x - fRes[1]*(c.p1.y - c.p0.y) - 2*fRes[0]*(c.p2.y - c.p1.y);
768 res.p2.x = c.p2.x - 2*fRes[1]*(c.p2.y - c.p1.y) - fRes[0]*(c.p3.y - c.p2.y);
769 res.p3.x = c.p3.x - 3*fRes[1]*(c.p3.y - c.p2.y);
771 res.p0.y = c.p0.y + 3*fRes[0]*(c.p1.x - c.p0.x);
772 res.p1.y = c.p1.y + 2*fRes[0]*(c.p2.x - c.p1.x) + fRes[1]*(c.p1.x - c.p0.x);
773 res.p2.y = c.p2.y + fRes[0]*(c.p3.x - c.p2.x) + 2*fRes[1]*(c.p2.x - c.p1.x);
774 res.p3.y = c.p3.y + 3*fRes[1]*(c.p3.x - c.p2.x);
777 // -----------------------------------------------------------------------------
779 bool Impl_calcSafeParams_focus( double& t1,
780 double& t2,
781 const Bezier& curve,
782 const Bezier& focus )
784 // now, we want to determine which normals of the original curve
785 // P(t) intersect with the focus curve F(t). The condition for
786 // this statement is P'(t)(P(t) - F) = 0, i.e. hodograph P'(t) and
787 // line through P(t) and F are perpendicular.
788 // If you expand this equation, you end up with something like
790 // (\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))
792 // Multiplying that out (as the scalar product is linear, we can
793 // extract some terms) yields:
795 // (P_i - F)^T n(P_{j+1} - P_j) B_i^n(t)B_j^{n-1}(t) + ...
797 // If we combine the B_i^n(t)B_j^{n-1}(t) product, we arrive at a
798 // Bernstein polynomial of degree 2n-1, as
800 // \binom{n}{i}(1-t)^{n-i}t^i) \binom{n-1}{j}(1-t)^{n-1-j}t^j) =
801 // \binom{n}{i}\binom{n-1}{j}(1-t)^{2n-1-i-j}t^{i+j}
803 // Thus, with the defining equation for a 2n-1 degree Bernstein
804 // polynomial
806 // \sum_{i=0}^{2n-1} d_i B_i^{2n-1}(t)
808 // the d_i are calculated as follows:
810 // 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)
813 // Okay, but F is now not a single point, but itself a curve
814 // F(u). Thus, for every value of u, we get a different 2n-1
815 // bezier curve from the above equation. Therefore, we have a
816 // tensor product bezier patch, with the following defining
817 // equation:
819 // d(t,u) = \sum_{i=0}^{2n-1} \sum_{j=0}^m B_i^{2n-1}(t) B_j^{m}(u) d_{ij}, where
820 // 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)
822 // as above, only that now F is one of the focus' control points.
824 // Note the difference in the binomial coefficients to the
825 // reference paper, these formulas most probably contained a typo.
827 // To determine, where D(t,u) is _not_ zero (these are the parts
828 // of the curve that don't share normals with the focus and can
829 // thus be safely clipped away), we project D(u,t) onto the
830 // (d(t,u), t) plane, determine the convex hull there and proceed
831 // as for the curve intersection part (projection is orthogonal to
832 // u axis, thus simply throw away u coordinate).
834 // \fallfac are so-called falling factorials (see Concrete
835 // Mathematics, p. 47 for a definition).
838 // now, for tensor product bezier curves, the convex hull property
839 // holds, too. Thus, we simply project the control points (t_{ij},
840 // u_{ij}, d_{ij}) onto the (t,d) plane and calculate the
841 // intersections of the convex hull with the t axis, as for the
842 // bezier clipping case.
845 // calc polygon of control points (t_{ij}, d_{ij}):
847 const int n( 3 ); // cubic bezier curves, as a matter of fact
848 const int i_card( 2*n );
849 const int j_card( n + 1 );
850 const int k_max( n-1 );
851 Polygon2D controlPolygon( i_card*j_card ); // vector of (t_{ij}, d_{ij}) in row-major order
853 int i, j, k, l; // variable notation from formulas above and Sederberg article
854 Point2D::value_type d;
855 for( i=0; i<i_card; ++i )
857 for( j=0; j<j_card; ++j )
859 // calc single d_{ij} sum:
860 for( d=0.0, k=::std::max(0,i-n); k<=k_max && k<=i; ++k )
862 l = i - k; // invariant: k + l = i
863 assert(k>=0 && k<=n-1); // k \in {0,...,n-1}
864 assert(l>=0 && l<=n); // l \in {0,...,n}
866 // TODO: find, document and assert proper limits for n and int's max_val.
867 // This becomes important should anybody wants to use
868 // this code for higher-than-cubic beziers
869 d += static_cast<double>(fallFac(n,l)*fallFac(n-1,k)*fac(i)) /
870 static_cast<double>(fac(l)*fac(k) * fallFac(2*n-1,i)) * n *
871 ( (curve[k+1].x - curve[k].x)*(curve[l].x - focus[j].x) + // dot product here
872 (curve[k+1].y - curve[k].y)*(curve[l].y - focus[j].y) );
875 // Note that the t_{ij} values are evenly spaced on the
876 // [0,1] interval, thus t_{ij}=i/(2n-1)
877 controlPolygon[ i*j_card + j ] = Point2D( i/(2.0*n-1.0), d );
881 #ifndef WITH_SAFEFOCUSPARAM_DETAILED_TEST
883 // calc safe parameter range, to determine [0,t1] and [t2,1] where
884 // no zero crossing is guaranteed.
885 return Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 );
887 #else
888 bool bRet( Impl_calcSafeParams( t1, t2, controlPolygon, 0.0, 0.0 ) );
890 Polygon2D convHull( convexHull( controlPolygon ) );
892 cout << "# convex hull testing (focus)" << endl
893 << "plot [t=0:1] ";
894 cout << "'-' using ($1):($2) title \"control polygon\" with lp, "
895 << "'-' using ($1):($2) title \"convex hull\" with lp" << endl;
897 unsigned int count;
898 for( count=0; count<controlPolygon.size(); ++count )
900 cout << controlPolygon[count].x << " " << controlPolygon[count].y << endl;
902 cout << controlPolygon[0].x << " " << controlPolygon[0].y << endl;
903 cout << "e" << endl;
905 for( count=0; count<convHull.size(); ++count )
907 cout << convHull[count].x << " " << convHull[count].y << endl;
909 cout << convHull[0].x << " " << convHull[0].y << endl;
910 cout << "e" << endl;
912 return bRet;
913 #endif
916 // -----------------------------------------------------------------------------
918 /** Calc all values t_i on c1, for which safeRanges functor does not
919 give a safe range on c1 and c2.
921 This method is the workhorse of the bezier clipping. Because c1
922 and c2 must be alternatingly tested against each other (first
923 determine safe parameter interval on c1 with regard to c2, then
924 the other way around), we call this method recursively with c1 and
925 c2 swapped.
927 @param result
928 Output iterator where the final t values are added to. If curves
929 don't intersect, nothing is added.
931 @param delta
932 Maximal allowed distance to true critical point (measured in the
933 original curve's coordinate system)
935 @param safeRangeFunctor
936 Functor object, that must provide the following operator():
937 bool safeRangeFunctor( double& t1,
938 double& t2,
939 const Bezier& c1_orig,
940 const Bezier& c1_part,
941 const Bezier& c2_orig,
942 const Bezier& c2_part );
943 This functor must calculate the safe ranges [0,t1] and [t2,1] on
944 c1_orig, where c1_orig is 'safe' from c2_part. If the whole
945 c1_orig is safe, false must be returned, true otherwise.
947 template <class Functor> void Impl_applySafeRanges_rec( ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > >& result,
948 double delta,
949 const Functor& safeRangeFunctor,
950 int recursionLevel,
951 const Bezier& c1_orig,
952 const Bezier& c1_part,
953 double last_t1_c1,
954 double last_t2_c1,
955 const Bezier& c2_orig,
956 const Bezier& c2_part,
957 double last_t1_c2,
958 double last_t2_c2 )
960 // check end condition
961 // ===================
963 // TODO: tidy up recursion handling. maybe put everything in a
964 // struct and swap that here at method entry
966 // TODO: Implement limit on recursion depth. Should that limit be
967 // reached, chances are that we're on a higher-order tangency. For
968 // this case, AW proposed to take the middle of the current
969 // interval, and to correct both curve's tangents at that new
970 // endpoint to be equal. That virtually generates a first-order
971 // tangency, and justifies to return a single intersection
972 // point. Otherwise, inside/outside test might fail here.
974 for( int i=0; i<recursionLevel; ++i ) cerr << " ";
975 if( recursionLevel % 2 )
977 cerr << "level: " << recursionLevel
978 << " t: "
979 << last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0
980 << ", c1: " << last_t1_c2 << " " << last_t2_c2
981 << ", c2: " << last_t1_c1 << " " << last_t2_c1
982 << endl;
984 else
986 cerr << "level: " << recursionLevel
987 << " t: "
988 << last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0
989 << ", c1: " << last_t1_c1 << " " << last_t2_c1
990 << ", c2: " << last_t1_c2 << " " << last_t2_c2
991 << endl;
994 // refine solution
995 // ===============
997 double t1_c1, t2_c1;
999 // Note: we first perform the clipping and only test for precision
1000 // sufficiency afterwards, since we want to exploit the fact that
1001 // Impl_calcClipRange returns false if the curves don't
1002 // intersect. We would have to check that separately for the end
1003 // condition, otherwise.
1005 // determine safe range on c1_orig
1006 if( safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part ) )
1008 // now, t1 and t2 are calculated on the original curve
1009 // (but against a fat line calculated from the subdivided
1010 // c2, namely c2_part). If the [t1,t2] range is outside
1011 // our current [last_t1,last_t2] range, we're done in this
1012 // branch - the curves no longer intersect.
1013 if( tolLessEqual(t1_c1, last_t2_c1) && tolGreaterEqual(t2_c1, last_t1_c1) )
1015 // As noted above, t1 and t2 are calculated on the
1016 // original curve, but against a fat line
1017 // calculated from the subdivided c2, namely
1018 // c2_part. Our domain to work on is
1019 // [last_t1,last_t2], on the other hand, so values
1020 // of [t1,t2] outside that range are irrelevant
1021 // here. Clip range appropriately.
1022 t1_c1 = ::std::max(t1_c1, last_t1_c1);
1023 t2_c1 = ::std::min(t2_c1, last_t2_c1);
1025 // TODO: respect delta
1026 // for now, end condition is just a fixed threshold on the t's
1028 // check end condition
1029 // ===================
1031 #if 1
1032 if( fabs(last_t2_c1 - last_t1_c1) < 0.0001 &&
1033 fabs(last_t2_c2 - last_t1_c2) < 0.0001 )
1034 #else
1035 if( fabs(last_t2_c1 - last_t1_c1) < 0.01 &&
1036 fabs(last_t2_c2 - last_t1_c2) < 0.01 )
1037 #endif
1039 // done. Add to result
1040 if( recursionLevel % 2 )
1042 // uneven level: have to swap the t's, since curves are swapped, too
1043 *result++ = ::std::make_pair( last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0,
1044 last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0 );
1046 else
1048 *result++ = ::std::make_pair( last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0,
1049 last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0 );
1052 #if 0
1053 //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, last_t1_c1, last_t2_c1 );
1054 printResultWithFinalCurves( c1_orig, c1_part, c2_orig, c2_part, t1_c1, t2_c1 );
1055 #else
1056 // calc focus curve of c2
1057 Bezier focus;
1058 Impl_calcFocus(focus, c2_part); // need to use subdivided c2
1060 safeRangeFunctor( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part );
1062 //printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, t1_c1, t2_c1 );
1063 printResultWithFinalCurves( c1_orig, c1_part, c2_orig, focus, last_t1_c1, last_t2_c1 );
1064 #endif
1066 else
1068 // heuristic: if parameter range is not reduced by at least
1069 // 20%, subdivide longest curve, and clip shortest against
1070 // both parts of longest
1071 // if( (last_t2_c1 - last_t1_c1 - t2_c1 + t1_c1) / (last_t2_c1 - last_t1_c1) < 0.2 )
1072 if( false )
1074 // subdivide and descend
1075 // =====================
1077 Bezier part1;
1078 Bezier part2;
1080 double intervalMiddle;
1082 if( last_t2_c1 - last_t1_c1 > last_t2_c2 - last_t1_c2 )
1084 // subdivide c1
1085 // ============
1087 intervalMiddle = last_t1_c1 + (last_t2_c1 - last_t1_c1)/2.0;
1089 // subdivide at the middle of the interval (as
1090 // we're not subdividing on the original
1091 // curve, this simply amounts to subdivision
1092 // at 0.5)
1093 Impl_deCasteljauAt( part1, part2, c1_part, 0.5 );
1095 // and descend recursively with swapped curves
1096 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1097 c2_orig, c2_part, last_t1_c2, last_t2_c2,
1098 c1_orig, part1, last_t1_c1, intervalMiddle );
1100 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1101 c2_orig, c2_part, last_t1_c2, last_t2_c2,
1102 c1_orig, part2, intervalMiddle, last_t2_c1 );
1104 else
1106 // subdivide c2
1107 // ============
1109 intervalMiddle = last_t1_c2 + (last_t2_c2 - last_t1_c2)/2.0;
1111 // subdivide at the middle of the interval (as
1112 // we're not subdividing on the original
1113 // curve, this simply amounts to subdivision
1114 // at 0.5)
1115 Impl_deCasteljauAt( part1, part2, c2_part, 0.5 );
1117 // and descend recursively with swapped curves
1118 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1119 c2_orig, part1, last_t1_c2, intervalMiddle,
1120 c1_orig, c1_part, last_t1_c1, last_t2_c1 );
1122 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1123 c2_orig, part2, intervalMiddle, last_t2_c2,
1124 c1_orig, c1_part, last_t1_c1, last_t2_c1 );
1127 else
1129 // apply calculated clip
1130 // =====================
1132 // clip safe ranges off c1_orig
1133 Bezier c1_part1;
1134 Bezier c1_part2;
1135 Bezier c1_part3;
1137 // subdivide at t1_c1
1138 Impl_deCasteljauAt( c1_part1, c1_part2, c1_orig, t1_c1 );
1140 // subdivide at t2_c1. As we're working on
1141 // c1_part2 now, we have to adapt t2_c1 since
1142 // we're no longer in the original parameter
1143 // interval. This is based on the following
1144 // assumption: t2_new = (t2-t1)/(1-t1), which
1145 // relates the t2 value into the new parameter
1146 // range [0,1] of c1_part2.
1147 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2_c1-t1_c1)/(1.0-t1_c1) );
1149 // descend with swapped curves and c1_part1 as the
1150 // remaining (middle) segment
1151 Impl_applySafeRanges_rec( result, delta, safeRangeFunctor, recursionLevel+1,
1152 c2_orig, c2_part, last_t1_c2, last_t2_c2,
1153 c1_orig, c1_part1, t1_c1, t2_c1 );
1160 // -----------------------------------------------------------------------------
1162 struct ClipBezierFunctor
1164 bool operator()( double& t1_c1,
1165 double& t2_c1,
1166 const Bezier& c1_orig,
1167 const Bezier& c1_part,
1168 const Bezier& c2_orig,
1169 const Bezier& c2_part ) const
1171 return Impl_calcClipRange( t1_c1, t2_c1, c1_orig, c1_part, c2_orig, c2_part );
1175 // -----------------------------------------------------------------------------
1177 struct BezierTangencyFunctor
1179 bool operator()( double& t1_c1,
1180 double& t2_c1,
1181 const Bezier& c1_orig,
1182 const Bezier& c1_part,
1183 const Bezier& c2_orig,
1184 const Bezier& c2_part ) const
1186 // calc focus curve of c2
1187 Bezier focus;
1188 Impl_calcFocus(focus, c2_part); // need to use subdivided c2
1189 // here, as the whole curve is
1190 // used for focus calculation
1192 // determine safe range on c1_orig
1193 bool bRet( Impl_calcSafeParams_focus( t1_c1, t2_c1,
1194 c1_orig, // use orig curve here, need t's on original curve
1195 focus ) );
1197 cerr << "range: " << t2_c1 - t1_c1 << ", ret: " << bRet << endl;
1199 return bRet;
1203 // -----------------------------------------------------------------------------
1205 /** Perform a bezier clip (curve against curve)
1207 @param result
1208 Output iterator where the final t values are added to. This
1209 iterator will remain empty, if there are no intersections.
1211 @param delta
1212 Maximal allowed distance to true intersection (measured in the
1213 original curve's coordinate system)
1215 void clipBezier( ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > >& result,
1216 double delta,
1217 const Bezier& c1,
1218 const Bezier& c2 )
1220 #if 0
1221 // first of all, determine list of collinear normals. Collinear
1222 // normals typically separate two intersections, thus, subdivide
1223 // at all collinear normal's t values beforehand. This will cater
1224 // for tangent intersections, where two or more intersections are
1225 // infinitesimally close together.
1227 // TODO: evaluate effects of higher-than-second-order
1228 // tangencies. Sederberg et al. state that collinear normal
1229 // algorithm then degrades quickly.
1231 ::std::vector< ::std::pair<double,double> > results;
1232 ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > > ii(results);
1234 Impl_calcCollinearNormals( ii, delta, 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
1236 // As Sederberg's collinear normal theorem is only sufficient, not
1237 // necessary for two intersections left and right, we've to test
1238 // all segments generated by the collinear normal algorithm
1239 // against each other. In other words, if the two curves are both
1240 // divided in a left and a right part, the collinear normal
1241 // theorem does _not_ state that the left part of curve 1 does not
1242 // e.g. intersect with the right part of curve 2.
1244 // divide c1 and c2 at collinear normal intersection points
1245 ::std::vector< Bezier > c1_segments( results.size()+1 );
1246 ::std::vector< Bezier > c2_segments( results.size()+1 );
1247 Bezier c1_remainder( c1 );
1248 Bezier c2_remainder( c2 );
1249 unsigned int i;
1250 for( i=0; i<results.size(); ++i )
1252 Bezier c1_part2;
1253 Impl_deCasteljauAt( c1_segments[i], c1_part2, c1_remainder, results[i].first );
1254 c1_remainder = c1_part2;
1256 Bezier c2_part2;
1257 Impl_deCasteljauAt( c2_segments[i], c2_part2, c2_remainder, results[i].second );
1258 c2_remainder = c2_part2;
1260 c1_segments[i] = c1_remainder;
1261 c2_segments[i] = c2_remainder;
1263 // now, c1/c2_segments contain all segments, then
1264 // clip every resulting segment against every other
1265 unsigned int c1_curr, c2_curr;
1266 for( c1_curr=0; c1_curr<c1_segments.size(); ++c1_curr )
1268 for( c2_curr=0; c2_curr<c2_segments.size(); ++c2_curr )
1270 if( c1_curr != c2_curr )
1272 Impl_clipBezier_rec(result, delta, 0,
1273 c1_segments[c1_curr], c1_segments[c1_curr],
1274 0.0, 1.0,
1275 c2_segments[c2_curr], c2_segments[c2_curr],
1276 0.0, 1.0);
1280 #else
1281 Impl_applySafeRanges_rec( result, delta, BezierTangencyFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
1282 //Impl_applySafeRanges_rec( result, delta, ClipBezierFunctor(), 0, c1, c1, 0.0, 1.0, c2, c2, 0.0, 1.0 );
1283 #endif
1284 // that's it, boys'n'girls!
1287 int main(int argc, const char *argv[])
1289 double curr_Offset( 0 );
1290 unsigned int i,j,k;
1292 Bezier someCurves[] =
1294 // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.0)},
1295 // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)},
1296 // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)}
1297 // {Point2D(0.25+1,0.5),Point2D(0.25+1,0.708333),Point2D(0.423611+1,0.916667),Point2D(0.770833+1,0.980324)},
1298 // {Point2D(0.0+1,0.0),Point2D(0.0+1,1.0),Point2D(1.0+1,1.0),Point2D(1.0+1,0.5)}
1300 // tangency1
1301 // {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)},
1302 // {Point2D(0.0,1.0),Point2D(0.1,1.0),Point2D(0.4,1.0),Point2D(0.5,1.0)}
1304 // {Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0),Point2D(1.0,0.5)},
1305 // {Point2D(0.60114,0.933091),Point2D(0.69461,0.969419),Point2D(0.80676,0.992976),Point2D(0.93756,0.998663)}
1306 // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)},
1307 // {Point2D(0.62712,0.828427),Point2D(0.76305,0.828507),Point2D(0.88555,0.77312),Point2D(0.95069,0.67325)}
1309 // clipping1
1310 // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
1311 // {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)}
1313 // tangency2
1314 // {Point2D(0.0,1.0),Point2D(3.5,1.0),Point2D(-2.5,0.0),Point2D(1.0,0.0)},
1315 // {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)}
1317 // tangency3
1318 // {Point2D(1.0,0.0),Point2D(0.0,0.0),Point2D(0.0,1.0),Point2D(1.0,1.0)},
1319 // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)}
1321 // tangency4
1322 // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)},
1323 // {Point2D(0.26,0.4),Point2D(0.25,0.5),Point2D(0.25,0.5),Point2D(0.26,0.6)}
1325 // tangency5
1326 // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
1327 // {Point2D(15.3621,0.00986464),Point2D(15.3683,0.0109389),Point2D(15.3682,0.0109315),Point2D(15.3621,0.00986464)}
1329 // tangency6
1330 // {Point2D(0.0,0.0),Point2D(0.0,3.5),Point2D(1.0,-2.5),Point2D(1.0,1.0)},
1331 // {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)}
1333 // tangency7
1334 // {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)},
1335 // {Point2D(15.3621,10.00986464),Point2D(15.3683,10.0109389),Point2D(15.3682,10.0109315),Point2D(15.3621,10.00986464)}
1337 // tangency Sederberg example
1338 {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)},
1339 {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)}
1341 // clipping2
1342 // {Point2D(-0.5,0.0),Point2D(0.5,0.0),Point2D(0.5,1.0),Point2D(-0.5,1.0)},
1343 // {Point2D(0.2575,0.4),Point2D(0.2475,0.5),Point2D(0.2475,0.5),Point2D(0.2575,0.6)}
1345 // {Point2D(0.0,0.1),Point2D(0.2,3.5),Point2D(1.0,-2.5),Point2D(1.1,1.2)},
1346 // {Point2D(0.0,1.0),Point2D(3.5,0.9),Point2D(-2.5,0.1),Point2D(1.1,0.2)}
1347 // {Point2D(0.0,0.1),Point2D(0.2,3.0),Point2D(1.0,-2.0),Point2D(1.1,1.2)},
1348 // {Point2D(0.627124+1,0.828427),Point2D(0.763048+1,0.828507),Point2D(0.885547+1,0.77312),Point2D(0.950692+1,0.67325)}
1349 // {Point2D(0.0,1.0),Point2D(3.0,0.9),Point2D(-2.0,0.1),Point2D(1.1,0.2)}
1350 // {Point2D(0.0,4.0),Point2D(0.1,5.0),Point2D(0.9,5.0),Point2D(1.0,4.0)},
1351 // {Point2D(0.0,0.0),Point2D(0.1,0.5),Point2D(0.9,0.5),Point2D(1.0,0.0)},
1352 // {Point2D(0.0,0.1),Point2D(0.1,1.5),Point2D(0.9,1.5),Point2D(1.0,0.1)},
1353 // {Point2D(0.0,-4.0),Point2D(0.1,-5.0),Point2D(0.9,-5.0),Point2D(1.0,-4.0)}
1356 // output gnuplot setup
1357 cout << "#!/usr/bin/gnuplot -persist" << endl
1358 << "#" << endl
1359 << "# automatically generated by bezierclip, don't change!" << endl
1360 << "#" << endl
1361 << "set parametric" << endl
1362 << "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
1363 << "bezd(p,q,r,s,t) = 3*(q-p)*(1-t)**2+6*(r-q)*(1-t)*t+3*(s-r)*t**2" << endl
1364 << "pointmarkx(c,t) = c-0.03*t" << endl
1365 << "pointmarky(c,t) = c+0.03*t" << endl
1366 << "linex(a,b,c,t) = a*-c + t*-b" << endl
1367 << "liney(a,b,c,t) = b*-c + t*a" << endl << endl
1368 << "# end of setup" << endl << endl;
1370 #ifdef WITH_CONVEXHULL_TEST
1371 // test convex hull algorithm
1372 const double convHull_xOffset( curr_Offset );
1373 curr_Offset += 20;
1374 cout << "# convex hull testing" << endl
1375 << "plot [t=0:1] ";
1376 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1378 Polygon2D aTestPoly(4);
1379 aTestPoly[0] = someCurves[i].p0;
1380 aTestPoly[1] = someCurves[i].p1;
1381 aTestPoly[2] = someCurves[i].p2;
1382 aTestPoly[3] = someCurves[i].p3;
1384 aTestPoly[0].x += convHull_xOffset;
1385 aTestPoly[1].x += convHull_xOffset;
1386 aTestPoly[2].x += convHull_xOffset;
1387 aTestPoly[3].x += convHull_xOffset;
1389 cout << " bez("
1390 << aTestPoly[0].x << ","
1391 << aTestPoly[1].x << ","
1392 << aTestPoly[2].x << ","
1393 << aTestPoly[3].x << ",t),bez("
1394 << aTestPoly[0].y << ","
1395 << aTestPoly[1].y << ","
1396 << aTestPoly[2].y << ","
1397 << aTestPoly[3].y << ",t), '-' using ($1):($2) title \"convex hull " << i << "\" with lp";
1399 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1400 cout << ",\\" << endl;
1401 else
1402 cout << endl;
1404 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1406 Polygon2D aTestPoly(4);
1407 aTestPoly[0] = someCurves[i].p0;
1408 aTestPoly[1] = someCurves[i].p1;
1409 aTestPoly[2] = someCurves[i].p2;
1410 aTestPoly[3] = someCurves[i].p3;
1412 aTestPoly[0].x += convHull_xOffset;
1413 aTestPoly[1].x += convHull_xOffset;
1414 aTestPoly[2].x += convHull_xOffset;
1415 aTestPoly[3].x += convHull_xOffset;
1417 Polygon2D convHull( convexHull(aTestPoly) );
1419 for( k=0; k<convHull.size(); ++k )
1421 cout << convHull[k].x << " " << convHull[k].y << endl;
1423 cout << convHull[0].x << " " << convHull[0].y << endl;
1424 cout << "e" << endl;
1426 #endif
1428 #ifdef WITH_MULTISUBDIVIDE_TEST
1429 // test convex hull algorithm
1430 const double multiSubdivide_xOffset( curr_Offset );
1431 curr_Offset += 20;
1432 cout << "# multi subdivide testing" << endl
1433 << "plot [t=0:1] ";
1434 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1436 Bezier c( someCurves[i] );
1437 Bezier c1_part1;
1438 Bezier c1_part2;
1439 Bezier c1_part3;
1441 c.p0.x += multiSubdivide_xOffset;
1442 c.p1.x += multiSubdivide_xOffset;
1443 c.p2.x += multiSubdivide_xOffset;
1444 c.p3.x += multiSubdivide_xOffset;
1446 const double t1( 0.1+i/(3.0*sizeof(someCurves)/sizeof(Bezier)) );
1447 const double t2( 0.9-i/(3.0*sizeof(someCurves)/sizeof(Bezier)) );
1449 // subdivide at t1
1450 Impl_deCasteljauAt( c1_part1, c1_part2, c, t1 );
1452 // subdivide at t2_c1. As we're working on
1453 // c1_part2 now, we have to adapt t2_c1 since
1454 // we're no longer in the original parameter
1455 // interval. This is based on the following
1456 // assumption: t2_new = (t2-t1)/(1-t1), which
1457 // relates the t2 value into the new parameter
1458 // range [0,1] of c1_part2.
1459 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
1461 // subdivide at t2
1462 Impl_deCasteljauAt( c1_part3, c1_part2, c, t2 );
1464 cout << " bez("
1465 << c1_part1.p0.x << ","
1466 << c1_part1.p1.x << ","
1467 << c1_part1.p2.x << ","
1468 << c1_part1.p3.x << ",t), bez("
1469 << c1_part1.p0.y+0.01 << ","
1470 << c1_part1.p1.y+0.01 << ","
1471 << c1_part1.p2.y+0.01 << ","
1472 << c1_part1.p3.y+0.01 << ",t) title \"middle " << i << "\", "
1473 << " bez("
1474 << c1_part2.p0.x << ","
1475 << c1_part2.p1.x << ","
1476 << c1_part2.p2.x << ","
1477 << c1_part2.p3.x << ",t), bez("
1478 << c1_part2.p0.y << ","
1479 << c1_part2.p1.y << ","
1480 << c1_part2.p2.y << ","
1481 << c1_part2.p3.y << ",t) title \"right " << i << "\", "
1482 << " bez("
1483 << c1_part3.p0.x << ","
1484 << c1_part3.p1.x << ","
1485 << c1_part3.p2.x << ","
1486 << c1_part3.p3.x << ",t), bez("
1487 << c1_part3.p0.y << ","
1488 << c1_part3.p1.y << ","
1489 << c1_part3.p2.y << ","
1490 << c1_part3.p3.y << ",t) title \"left " << i << "\"";
1493 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1494 cout << ",\\" << endl;
1495 else
1496 cout << endl;
1498 #endif
1500 #ifdef WITH_FATLINE_TEST
1501 // test fatline algorithm
1502 const double fatLine_xOffset( curr_Offset );
1503 curr_Offset += 20;
1504 cout << "# fat line testing" << endl
1505 << "plot [t=0:1] ";
1506 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1508 Bezier c( someCurves[i] );
1510 c.p0.x += fatLine_xOffset;
1511 c.p1.x += fatLine_xOffset;
1512 c.p2.x += fatLine_xOffset;
1513 c.p3.x += fatLine_xOffset;
1515 FatLine line;
1517 Impl_calcFatLine(line, c);
1519 cout << " bez("
1520 << c.p0.x << ","
1521 << c.p1.x << ","
1522 << c.p2.x << ","
1523 << c.p3.x << ",t), bez("
1524 << c.p0.y << ","
1525 << c.p1.y << ","
1526 << c.p2.y << ","
1527 << c.p3.y << ",t) title \"bezier " << i << "\", linex("
1528 << line.a << ","
1529 << line.b << ","
1530 << line.c << ",t), liney("
1531 << line.a << ","
1532 << line.b << ","
1533 << line.c << ",t) title \"fat line (center) on " << i << "\", linex("
1534 << line.a << ","
1535 << line.b << ","
1536 << line.c-line.dMin << ",t), liney("
1537 << line.a << ","
1538 << line.b << ","
1539 << line.c-line.dMin << ",t) title \"fat line (min) on " << i << "\", linex("
1540 << line.a << ","
1541 << line.b << ","
1542 << line.c-line.dMax << ",t), liney("
1543 << line.a << ","
1544 << line.b << ","
1545 << line.c-line.dMax << ",t) title \"fat line (max) on " << i << "\"";
1547 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1548 cout << ",\\" << endl;
1549 else
1550 cout << endl;
1552 #endif
1554 #ifdef WITH_CALCFOCUS_TEST
1555 // test focus curve algorithm
1556 const double focus_xOffset( curr_Offset );
1557 curr_Offset += 20;
1558 cout << "# focus line testing" << endl
1559 << "plot [t=0:1] ";
1560 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1562 Bezier c( someCurves[i] );
1564 c.p0.x += focus_xOffset;
1565 c.p1.x += focus_xOffset;
1566 c.p2.x += focus_xOffset;
1567 c.p3.x += focus_xOffset;
1569 // calc focus curve
1570 Bezier focus;
1571 Impl_calcFocus(focus, c);
1573 cout << " bez("
1574 << c.p0.x << ","
1575 << c.p1.x << ","
1576 << c.p2.x << ","
1577 << c.p3.x << ",t), bez("
1578 << c.p0.y << ","
1579 << c.p1.y << ","
1580 << c.p2.y << ","
1581 << c.p3.y << ",t) title \"bezier " << i << "\", bez("
1582 << focus.p0.x << ","
1583 << focus.p1.x << ","
1584 << focus.p2.x << ","
1585 << focus.p3.x << ",t), bez("
1586 << focus.p0.y << ","
1587 << focus.p1.y << ","
1588 << focus.p2.y << ","
1589 << focus.p3.y << ",t) title \"focus " << i << "\"";
1592 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1593 cout << ",\\" << endl;
1594 else
1595 cout << endl;
1597 #endif
1599 #ifdef WITH_SAFEPARAMBASE_TEST
1600 // test safe params base method
1601 double safeParamsBase_xOffset( curr_Offset );
1602 cout << "# safe param base method testing" << endl
1603 << "plot [t=0:1] ";
1604 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1606 Bezier c( someCurves[i] );
1608 c.p0.x += safeParamsBase_xOffset;
1609 c.p1.x += safeParamsBase_xOffset;
1610 c.p2.x += safeParamsBase_xOffset;
1611 c.p3.x += safeParamsBase_xOffset;
1613 Polygon2D poly(4);
1614 poly[0] = c.p0;
1615 poly[1] = c.p1;
1616 poly[2] = c.p2;
1617 poly[3] = c.p3;
1619 double t1, t2;
1621 bool bRet( Impl_calcSafeParams( t1, t2, poly, 0, 1 ) );
1623 Polygon2D convHull( convexHull( poly ) );
1625 cout << " bez("
1626 << poly[0].x << ","
1627 << poly[1].x << ","
1628 << poly[2].x << ","
1629 << poly[3].x << ",t),bez("
1630 << poly[0].y << ","
1631 << poly[1].y << ","
1632 << poly[2].y << ","
1633 << poly[3].y << ",t), "
1634 << "t+" << safeParamsBase_xOffset << ", 0, "
1635 << "t+" << safeParamsBase_xOffset << ", 1, ";
1636 if( bRet )
1638 cout << t1+safeParamsBase_xOffset << ", t, "
1639 << t2+safeParamsBase_xOffset << ", t, ";
1641 cout << "'-' using ($1):($2) title \"control polygon\" with lp, "
1642 << "'-' using ($1):($2) title \"convex hull\" with lp";
1644 if( i+1<sizeof(someCurves)/sizeof(Bezier) )
1645 cout << ",\\" << endl;
1646 else
1647 cout << endl;
1649 safeParamsBase_xOffset += 2;
1652 safeParamsBase_xOffset = curr_Offset;
1653 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1655 Bezier c( someCurves[i] );
1657 c.p0.x += safeParamsBase_xOffset;
1658 c.p1.x += safeParamsBase_xOffset;
1659 c.p2.x += safeParamsBase_xOffset;
1660 c.p3.x += safeParamsBase_xOffset;
1662 Polygon2D poly(4);
1663 poly[0] = c.p0;
1664 poly[1] = c.p1;
1665 poly[2] = c.p2;
1666 poly[3] = c.p3;
1668 double t1, t2;
1670 Impl_calcSafeParams( t1, t2, poly, 0, 1 );
1672 Polygon2D convHull( convexHull( poly ) );
1674 unsigned int k;
1675 for( k=0; k<poly.size(); ++k )
1677 cout << poly[k].x << " " << poly[k].y << endl;
1679 cout << poly[0].x << " " << poly[0].y << endl;
1680 cout << "e" << endl;
1682 for( k=0; k<convHull.size(); ++k )
1684 cout << convHull[k].x << " " << convHull[k].y << endl;
1686 cout << convHull[0].x << " " << convHull[0].y << endl;
1687 cout << "e" << endl;
1689 safeParamsBase_xOffset += 2;
1691 curr_Offset += 20;
1692 #endif
1694 #ifdef WITH_SAFEPARAMS_TEST
1695 // test safe parameter range algorithm
1696 const double safeParams_xOffset( curr_Offset );
1697 curr_Offset += 20;
1698 cout << "# safe param range testing" << endl
1699 << "plot [t=0.0:1.0] ";
1700 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1702 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
1704 Bezier c1( someCurves[i] );
1705 Bezier c2( someCurves[j] );
1707 c1.p0.x += safeParams_xOffset;
1708 c1.p1.x += safeParams_xOffset;
1709 c1.p2.x += safeParams_xOffset;
1710 c1.p3.x += safeParams_xOffset;
1711 c2.p0.x += safeParams_xOffset;
1712 c2.p1.x += safeParams_xOffset;
1713 c2.p2.x += safeParams_xOffset;
1714 c2.p3.x += safeParams_xOffset;
1716 double t1, t2;
1718 if( Impl_calcClipRange(t1, t2, c1, c1, c2, c2) )
1720 // clip safe ranges off c1
1721 Bezier c1_part1;
1722 Bezier c1_part2;
1723 Bezier c1_part3;
1725 // subdivide at t1_c1
1726 Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 );
1727 // subdivide at t2_c1
1728 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
1730 // output remaining segment (c1_part1)
1732 cout << " bez("
1733 << c1.p0.x << ","
1734 << c1.p1.x << ","
1735 << c1.p2.x << ","
1736 << c1.p3.x << ",t),bez("
1737 << c1.p0.y << ","
1738 << c1.p1.y << ","
1739 << c1.p2.y << ","
1740 << c1.p3.y << ",t), bez("
1741 << c2.p0.x << ","
1742 << c2.p1.x << ","
1743 << c2.p2.x << ","
1744 << c2.p3.x << ",t),bez("
1745 << c2.p0.y << ","
1746 << c2.p1.y << ","
1747 << c2.p2.y << ","
1748 << c2.p3.y << ",t), bez("
1749 << c1_part1.p0.x << ","
1750 << c1_part1.p1.x << ","
1751 << c1_part1.p2.x << ","
1752 << c1_part1.p3.x << ",t),bez("
1753 << c1_part1.p0.y << ","
1754 << c1_part1.p1.y << ","
1755 << c1_part1.p2.y << ","
1756 << c1_part1.p3.y << ",t)";
1758 if( i+2<sizeof(someCurves)/sizeof(Bezier) )
1759 cout << ",\\" << endl;
1760 else
1761 cout << endl;
1765 #endif
1767 #ifdef WITH_SAFEPARAM_DETAILED_TEST
1768 // test safe parameter range algorithm
1769 const double safeParams2_xOffset( curr_Offset );
1770 curr_Offset += 20;
1771 if( sizeof(someCurves)/sizeof(Bezier) > 1 )
1773 Bezier c1( someCurves[0] );
1774 Bezier c2( someCurves[1] );
1776 c1.p0.x += safeParams2_xOffset;
1777 c1.p1.x += safeParams2_xOffset;
1778 c1.p2.x += safeParams2_xOffset;
1779 c1.p3.x += safeParams2_xOffset;
1780 c2.p0.x += safeParams2_xOffset;
1781 c2.p1.x += safeParams2_xOffset;
1782 c2.p2.x += safeParams2_xOffset;
1783 c2.p3.x += safeParams2_xOffset;
1785 double t1, t2;
1787 // output happens here
1788 Impl_calcClipRange(t1, t2, c1, c1, c2, c2);
1790 #endif
1792 #ifdef WITH_SAFEFOCUSPARAM_TEST
1793 // test safe parameter range from focus algorithm
1794 const double safeParamsFocus_xOffset( curr_Offset );
1795 curr_Offset += 20;
1796 cout << "# safe param range from focus testing" << endl
1797 << "plot [t=0.0:1.0] ";
1798 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1800 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
1802 Bezier c1( someCurves[i] );
1803 Bezier c2( someCurves[j] );
1805 c1.p0.x += safeParamsFocus_xOffset;
1806 c1.p1.x += safeParamsFocus_xOffset;
1807 c1.p2.x += safeParamsFocus_xOffset;
1808 c1.p3.x += safeParamsFocus_xOffset;
1809 c2.p0.x += safeParamsFocus_xOffset;
1810 c2.p1.x += safeParamsFocus_xOffset;
1811 c2.p2.x += safeParamsFocus_xOffset;
1812 c2.p3.x += safeParamsFocus_xOffset;
1814 double t1, t2;
1816 Bezier focus;
1817 #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
1818 #if 0
1820 // clip safe ranges off c1_orig
1821 Bezier c1_part1;
1822 Bezier c1_part2;
1823 Bezier c1_part3;
1825 // subdivide at t1_c1
1826 Impl_deCasteljauAt( c1_part1, c1_part2, c2, 0.30204 );
1828 // subdivide at t2_c1. As we're working on
1829 // c1_part2 now, we have to adapt t2_c1 since
1830 // we're no longer in the original parameter
1831 // interval. This is based on the following
1832 // assumption: t2_new = (t2-t1)/(1-t1), which
1833 // relates the t2 value into the new parameter
1834 // range [0,1] of c1_part2.
1835 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (0.57151-0.30204)/(1.0-0.30204) );
1837 c2 = c1_part1;
1838 Impl_calcFocus( focus, c2 );
1840 #else
1841 Impl_calcFocus( focus, c2 );
1842 #endif
1843 #else
1844 focus = c2;
1845 #endif
1846 // determine safe range on c1
1847 bool bRet( Impl_calcSafeParams_focus( t1, t2,
1848 c1, focus ) );
1850 cerr << "t1: " << t1 << ", t2: " << t2 << endl;
1852 // clip safe ranges off c1
1853 Bezier c1_part1;
1854 Bezier c1_part2;
1855 Bezier c1_part3;
1857 // subdivide at t1_c1
1858 Impl_deCasteljauAt( c1_part1, c1_part2, c1, t1 );
1859 // subdivide at t2_c1
1860 Impl_deCasteljauAt( c1_part1, c1_part3, c1_part2, (t2-t1)/(1.0-t1) );
1862 // output remaining segment (c1_part1)
1864 cout << " bez("
1865 << c1.p0.x << ","
1866 << c1.p1.x << ","
1867 << c1.p2.x << ","
1868 << c1.p3.x << ",t),bez("
1869 << c1.p0.y << ","
1870 << c1.p1.y << ","
1871 << c1.p2.y << ","
1872 << c1.p3.y << ",t) title \"c1\", "
1873 #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
1874 << "bez("
1875 << c2.p0.x << ","
1876 << c2.p1.x << ","
1877 << c2.p2.x << ","
1878 << c2.p3.x << ",t),bez("
1879 << c2.p0.y << ","
1880 << c2.p1.y << ","
1881 << c2.p2.y << ","
1882 << c2.p3.y << ",t) title \"c2\", "
1883 << "bez("
1884 << focus.p0.x << ","
1885 << focus.p1.x << ","
1886 << focus.p2.x << ","
1887 << focus.p3.x << ",t),bez("
1888 << focus.p0.y << ","
1889 << focus.p1.y << ","
1890 << focus.p2.y << ","
1891 << focus.p3.y << ",t) title \"focus\"";
1892 #else
1893 << "bez("
1894 << c2.p0.x << ","
1895 << c2.p1.x << ","
1896 << c2.p2.x << ","
1897 << c2.p3.x << ",t),bez("
1898 << c2.p0.y << ","
1899 << c2.p1.y << ","
1900 << c2.p2.y << ","
1901 << c2.p3.y << ",t) title \"focus\"";
1902 #endif
1903 if( bRet )
1905 cout << ", bez("
1906 << c1_part1.p0.x << ","
1907 << c1_part1.p1.x << ","
1908 << c1_part1.p2.x << ","
1909 << c1_part1.p3.x << ",t),bez("
1910 << c1_part1.p0.y+0.01 << ","
1911 << c1_part1.p1.y+0.01 << ","
1912 << c1_part1.p2.y+0.01 << ","
1913 << c1_part1.p3.y+0.01 << ",t) title \"part\"";
1916 if( i+2<sizeof(someCurves)/sizeof(Bezier) )
1917 cout << ",\\" << endl;
1918 else
1919 cout << endl;
1922 #endif
1924 #ifdef WITH_SAFEFOCUSPARAM_DETAILED_TEST
1925 // test safe parameter range algorithm
1926 const double safeParams3_xOffset( curr_Offset );
1927 curr_Offset += 20;
1928 if( sizeof(someCurves)/sizeof(Bezier) > 1 )
1930 Bezier c1( someCurves[0] );
1931 Bezier c2( someCurves[1] );
1933 c1.p0.x += safeParams3_xOffset;
1934 c1.p1.x += safeParams3_xOffset;
1935 c1.p2.x += safeParams3_xOffset;
1936 c1.p3.x += safeParams3_xOffset;
1937 c2.p0.x += safeParams3_xOffset;
1938 c2.p1.x += safeParams3_xOffset;
1939 c2.p2.x += safeParams3_xOffset;
1940 c2.p3.x += safeParams3_xOffset;
1942 double t1, t2;
1944 Bezier focus;
1945 #ifdef WITH_SAFEFOCUSPARAM_CALCFOCUS
1946 Impl_calcFocus( focus, c2 );
1947 #else
1948 focus = c2;
1949 #endif
1951 // determine safe range on c1, output happens here
1952 Impl_calcSafeParams_focus( t1, t2,
1953 c1, focus );
1955 #endif
1957 #ifdef WITH_BEZIERCLIP_TEST
1958 ::std::vector< ::std::pair<double, double> > result;
1959 ::std::back_insert_iterator< ::std::vector< ::std::pair<double, double> > > ii(result);
1961 // test full bezier clipping
1962 const double bezierClip_xOffset( curr_Offset );
1963 curr_Offset += 20;
1964 cout << endl << endl << "# bezier clip testing" << endl
1965 << "plot [t=0:1] ";
1966 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
1968 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
1970 Bezier c1( someCurves[i] );
1971 Bezier c2( someCurves[j] );
1973 c1.p0.x += bezierClip_xOffset;
1974 c1.p1.x += bezierClip_xOffset;
1975 c1.p2.x += bezierClip_xOffset;
1976 c1.p3.x += bezierClip_xOffset;
1977 c2.p0.x += bezierClip_xOffset;
1978 c2.p1.x += bezierClip_xOffset;
1979 c2.p2.x += bezierClip_xOffset;
1980 c2.p3.x += bezierClip_xOffset;
1982 cout << " bez("
1983 << c1.p0.x << ","
1984 << c1.p1.x << ","
1985 << c1.p2.x << ","
1986 << c1.p3.x << ",t),bez("
1987 << c1.p0.y << ","
1988 << c1.p1.y << ","
1989 << c1.p2.y << ","
1990 << c1.p3.y << ",t), bez("
1991 << c2.p0.x << ","
1992 << c2.p1.x << ","
1993 << c2.p2.x << ","
1994 << c2.p3.x << ",t),bez("
1995 << c2.p0.y << ","
1996 << c2.p1.y << ","
1997 << c2.p2.y << ","
1998 << c2.p3.y << ",t), '-' using (bez("
1999 << c1.p0.x << ","
2000 << c1.p1.x << ","
2001 << c1.p2.x << ","
2002 << c1.p3.x
2003 << ",$1)):(bez("
2004 << c1.p0.y << ","
2005 << c1.p1.y << ","
2006 << c1.p2.y << ","
2007 << c1.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << i << ")\", "
2008 << " '-' using (bez("
2009 << c2.p0.x << ","
2010 << c2.p1.x << ","
2011 << c2.p2.x << ","
2012 << c2.p3.x
2013 << ",$1)):(bez("
2014 << c2.p0.y << ","
2015 << c2.p1.y << ","
2016 << c2.p2.y << ","
2017 << c2.p3.y << ",$1)) title \"bezier " << i << " clipped against " << j << " (t on " << j << ")\"";
2019 if( i+2<sizeof(someCurves)/sizeof(Bezier) )
2020 cout << ",\\" << endl;
2021 else
2022 cout << endl;
2025 for( i=0; i<sizeof(someCurves)/sizeof(Bezier); ++i )
2027 for( j=i+1; j<sizeof(someCurves)/sizeof(Bezier); ++j )
2029 result.clear();
2030 Bezier c1( someCurves[i] );
2031 Bezier c2( someCurves[j] );
2033 c1.p0.x += bezierClip_xOffset;
2034 c1.p1.x += bezierClip_xOffset;
2035 c1.p2.x += bezierClip_xOffset;
2036 c1.p3.x += bezierClip_xOffset;
2037 c2.p0.x += bezierClip_xOffset;
2038 c2.p1.x += bezierClip_xOffset;
2039 c2.p2.x += bezierClip_xOffset;
2040 c2.p3.x += bezierClip_xOffset;
2042 clipBezier( ii, 0.00001, c1, c2 );
2044 for( k=0; k<result.size(); ++k )
2046 cout << result[k].first << endl;
2048 cout << "e" << endl;
2050 for( k=0; k<result.size(); ++k )
2052 cout << result[k].second << endl;
2054 cout << "e" << endl;
2057 #endif
2059 return 0;