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