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