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