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