1//---------------------------------------------------------------------------- 2// Anti-Grain Geometry - Version 2.4 3// Copyright (C) 2002-2005 Maxim Shemanarev (http://www.antigrain.com) 4// 5// Permission to copy, use, modify, sell and distribute this software 6// is granted provided this copyright notice appears in all copies. 7// This software is provided "as is" without express or implied 8// warranty, and with no claim as to its suitability for any purpose. 9// 10//---------------------------------------------------------------------------- 11// Contact: mcseem@antigrain.com 12// mcseemagg@yahoo.com 13// http://www.antigrain.com 14//---------------------------------------------------------------------------- 15 16#include <math.h> 17#include "agg_curves.h" 18#include "agg_math.h" 19 20namespace agg 21{ 22 23 //------------------------------------------------------------------------ 24 const double curve_distance_epsilon = 1e-30; 25 const double curve_collinearity_epsilon = 1e-30; 26 const double curve_angle_tolerance_epsilon = 0.01; 27 enum curve_recursion_limit_e { curve_recursion_limit = 32 }; 28 29 30 31 //------------------------------------------------------------------------ 32 void curve3_inc::approximation_scale(double s) 33 { 34 m_scale = s; 35 } 36 37 //------------------------------------------------------------------------ 38 double curve3_inc::approximation_scale() const 39 { 40 return m_scale; 41 } 42 43 //------------------------------------------------------------------------ 44 void curve3_inc::init(double x1, double y1, 45 double x2, double y2, 46 double x3, double y3) 47 { 48 m_start_x = x1; 49 m_start_y = y1; 50 m_end_x = x3; 51 m_end_y = y3; 52 53 double dx1 = x2 - x1; 54 double dy1 = y2 - y1; 55 double dx2 = x3 - x2; 56 double dy2 = y3 - y2; 57 58 double len = sqrt(dx1 * dx1 + dy1 * dy1) + sqrt(dx2 * dx2 + dy2 * dy2); 59 60 m_num_steps = uround(len * 0.25 * m_scale); 61 62 if(m_num_steps < 4) 63 { 64 m_num_steps = 4; 65 } 66 67 double subdivide_step = 1.0 / m_num_steps; 68 double subdivide_step2 = subdivide_step * subdivide_step; 69 70 double tmpx = (x1 - x2 * 2.0 + x3) * subdivide_step2; 71 double tmpy = (y1 - y2 * 2.0 + y3) * subdivide_step2; 72 73 m_saved_fx = m_fx = x1; 74 m_saved_fy = m_fy = y1; 75 76 m_saved_dfx = m_dfx = tmpx + (x2 - x1) * (2.0 * subdivide_step); 77 m_saved_dfy = m_dfy = tmpy + (y2 - y1) * (2.0 * subdivide_step); 78 79 m_ddfx = tmpx * 2.0; 80 m_ddfy = tmpy * 2.0; 81 82 m_step = m_num_steps; 83 } 84 85 //------------------------------------------------------------------------ 86 void curve3_inc::rewind(unsigned) 87 { 88 if(m_num_steps == 0) 89 { 90 m_step = -1; 91 return; 92 } 93 m_step = m_num_steps; 94 m_fx = m_saved_fx; 95 m_fy = m_saved_fy; 96 m_dfx = m_saved_dfx; 97 m_dfy = m_saved_dfy; 98 } 99 100 //------------------------------------------------------------------------ 101 unsigned curve3_inc::vertex(double* x, double* y) 102 { 103 if(m_step < 0) return path_cmd_stop; 104 if(m_step == m_num_steps) 105 { 106 *x = m_start_x; 107 *y = m_start_y; 108 --m_step; 109 return path_cmd_move_to; 110 } 111 if(m_step == 0) 112 { 113 *x = m_end_x; 114 *y = m_end_y; 115 --m_step; 116 return path_cmd_line_to; 117 } 118 m_fx += m_dfx; 119 m_fy += m_dfy; 120 m_dfx += m_ddfx; 121 m_dfy += m_ddfy; 122 *x = m_fx; 123 *y = m_fy; 124 --m_step; 125 return path_cmd_line_to; 126 } 127 128 //------------------------------------------------------------------------ 129 void curve3_div::init(double x1, double y1, 130 double x2, double y2, 131 double x3, double y3) 132 { 133 m_points.remove_all(); 134 m_distance_tolerance_square = 0.5 / m_approximation_scale; 135 m_distance_tolerance_square *= m_distance_tolerance_square; 136 bezier(x1, y1, x2, y2, x3, y3); 137 m_count = 0; 138 } 139 140 //------------------------------------------------------------------------ 141 void curve3_div::recursive_bezier(double x1, double y1, 142 double x2, double y2, 143 double x3, double y3, 144 unsigned level) 145 { 146 if(level > curve_recursion_limit) 147 { 148 return; 149 } 150 151 // Calculate all the mid-points of the line segments 152 //---------------------- 153 double x12 = (x1 + x2) / 2; 154 double y12 = (y1 + y2) / 2; 155 double x23 = (x2 + x3) / 2; 156 double y23 = (y2 + y3) / 2; 157 double x123 = (x12 + x23) / 2; 158 double y123 = (y12 + y23) / 2; 159 160 double dx = x3-x1; 161 double dy = y3-y1; 162 double d = fabs(((x2 - x3) * dy - (y2 - y3) * dx)); 163 double da; 164 165 if(d > curve_collinearity_epsilon) 166 { 167 // Regular case 168 //----------------- 169 if(d * d <= m_distance_tolerance_square * (dx*dx + dy*dy)) 170 { 171 // If the curvature doesn't exceed the distance_tolerance value 172 // we tend to finish subdivisions. 173 //---------------------- 174 if(m_angle_tolerance < curve_angle_tolerance_epsilon) 175 { 176 m_points.add(point_d(x123, y123)); 177 return; 178 } 179 180 // Angle & Cusp Condition 181 //---------------------- 182 da = fabs(atan2(y3 - y2, x3 - x2) - atan2(y2 - y1, x2 - x1)); 183 if(da >= pi) da = 2*pi - da; 184 185 if(da < m_angle_tolerance) 186 { 187 // Finally we can stop the recursion 188 //---------------------- 189 m_points.add(point_d(x123, y123)); 190 return; 191 } 192 } 193 } 194 else 195 { 196 // Collinear case 197 //------------------ 198 da = dx*dx + dy*dy; 199 if(da == 0) 200 { 201 d = calc_sq_distance(x1, y1, x2, y2); 202 } 203 else 204 { 205 d = ((x2 - x1)*dx + (y2 - y1)*dy) / da; 206 if(d > 0 && d < 1) 207 { 208 // Simple collinear case, 1---2---3 209 // We can leave just two endpoints 210 return; 211 } 212 if(d <= 0) d = calc_sq_distance(x2, y2, x1, y1); 213 else if(d >= 1) d = calc_sq_distance(x2, y2, x3, y3); 214 else d = calc_sq_distance(x2, y2, x1 + d*dx, y1 + d*dy); 215 } 216 if(d < m_distance_tolerance_square) 217 { 218 m_points.add(point_d(x2, y2)); 219 return; 220 } 221 } 222 223 // Continue subdivision 224 //---------------------- 225 recursive_bezier(x1, y1, x12, y12, x123, y123, level + 1); 226 recursive_bezier(x123, y123, x23, y23, x3, y3, level + 1); 227 } 228 229 //------------------------------------------------------------------------ 230 void curve3_div::bezier(double x1, double y1, 231 double x2, double y2, 232 double x3, double y3) 233 { 234 m_points.add(point_d(x1, y1)); 235 recursive_bezier(x1, y1, x2, y2, x3, y3, 0); 236 m_points.add(point_d(x3, y3)); 237 } 238 239 240 241 242 243 //------------------------------------------------------------------------ 244 void curve4_inc::approximation_scale(double s) 245 { 246 m_scale = s; 247 } 248 249 //------------------------------------------------------------------------ 250 double curve4_inc::approximation_scale() const 251 { 252 return m_scale; 253 } 254 255 //------------------------------------------------------------------------ 256 static double MSC60_fix_ICE(double v) { return v; } 257 258 //------------------------------------------------------------------------ 259 void curve4_inc::init(double x1, double y1, 260 double x2, double y2, 261 double x3, double y3, 262 double x4, double y4) 263 { 264 m_start_x = x1; 265 m_start_y = y1; 266 m_end_x = x4; 267 m_end_y = y4; 268 269 double dx1 = x2 - x1; 270 double dy1 = y2 - y1; 271 double dx2 = x3 - x2; 272 double dy2 = y3 - y2; 273 double dx3 = x4 - x3; 274 double dy3 = y4 - y3; 275 276 double len = (sqrt(dx1 * dx1 + dy1 * dy1) + 277 sqrt(dx2 * dx2 + dy2 * dy2) + 278 sqrt(dx3 * dx3 + dy3 * dy3)) * 0.25 * m_scale; 279 280#if defined(_MSC_VER) && _MSC_VER <= 1200 281 m_num_steps = uround(MSC60_fix_ICE(len)); 282#else 283 m_num_steps = uround(len); 284#endif 285 286 if(m_num_steps < 4) 287 { 288 m_num_steps = 4; 289 } 290 291 double subdivide_step = 1.0 / m_num_steps; 292 double subdivide_step2 = subdivide_step * subdivide_step; 293 double subdivide_step3 = subdivide_step * subdivide_step * subdivide_step; 294 295 double pre1 = 3.0 * subdivide_step; 296 double pre2 = 3.0 * subdivide_step2; 297 double pre4 = 6.0 * subdivide_step2; 298 double pre5 = 6.0 * subdivide_step3; 299 300 double tmp1x = x1 - x2 * 2.0 + x3; 301 double tmp1y = y1 - y2 * 2.0 + y3; 302 303 double tmp2x = (x2 - x3) * 3.0 - x1 + x4; 304 double tmp2y = (y2 - y3) * 3.0 - y1 + y4; 305 306 m_saved_fx = m_fx = x1; 307 m_saved_fy = m_fy = y1; 308 309 m_saved_dfx = m_dfx = (x2 - x1) * pre1 + tmp1x * pre2 + tmp2x * subdivide_step3; 310 m_saved_dfy = m_dfy = (y2 - y1) * pre1 + tmp1y * pre2 + tmp2y * subdivide_step3; 311 312 m_saved_ddfx = m_ddfx = tmp1x * pre4 + tmp2x * pre5; 313 m_saved_ddfy = m_ddfy = tmp1y * pre4 + tmp2y * pre5; 314 315 m_dddfx = tmp2x * pre5; 316 m_dddfy = tmp2y * pre5; 317 318 m_step = m_num_steps; 319 } 320 321 //------------------------------------------------------------------------ 322 void curve4_inc::rewind(unsigned) 323 { 324 if(m_num_steps == 0) 325 { 326 m_step = -1; 327 return; 328 } 329 m_step = m_num_steps; 330 m_fx = m_saved_fx; 331 m_fy = m_saved_fy; 332 m_dfx = m_saved_dfx; 333 m_dfy = m_saved_dfy; 334 m_ddfx = m_saved_ddfx; 335 m_ddfy = m_saved_ddfy; 336 } 337 338 //------------------------------------------------------------------------ 339 unsigned curve4_inc::vertex(double* x, double* y) 340 { 341 if(m_step < 0) return path_cmd_stop; 342 if(m_step == m_num_steps) 343 { 344 *x = m_start_x; 345 *y = m_start_y; 346 --m_step; 347 return path_cmd_move_to; 348 } 349 350 if(m_step == 0) 351 { 352 *x = m_end_x; 353 *y = m_end_y; 354 --m_step; 355 return path_cmd_line_to; 356 } 357 358 m_fx += m_dfx; 359 m_fy += m_dfy; 360 m_dfx += m_ddfx; 361 m_dfy += m_ddfy; 362 m_ddfx += m_dddfx; 363 m_ddfy += m_dddfy; 364 365 *x = m_fx; 366 *y = m_fy; 367 --m_step; 368 return path_cmd_line_to; 369 } 370 371 372 373 374 //------------------------------------------------------------------------ 375 void curve4_div::init(double x1, double y1, 376 double x2, double y2, 377 double x3, double y3, 378 double x4, double y4) 379 { 380 m_points.remove_all(); 381 m_distance_tolerance_square = 0.5 / m_approximation_scale; 382 m_distance_tolerance_square *= m_distance_tolerance_square; 383 bezier(x1, y1, x2, y2, x3, y3, x4, y4); 384 m_count = 0; 385 } 386 387 //------------------------------------------------------------------------ 388 void curve4_div::recursive_bezier(double x1, double y1, 389 double x2, double y2, 390 double x3, double y3, 391 double x4, double y4, 392 unsigned level) 393 { 394 if(level > curve_recursion_limit) 395 { 396 return; 397 } 398 399 // Calculate all the mid-points of the line segments 400 //---------------------- 401 double x12 = (x1 + x2) / 2; 402 double y12 = (y1 + y2) / 2; 403 double x23 = (x2 + x3) / 2; 404 double y23 = (y2 + y3) / 2; 405 double x34 = (x3 + x4) / 2; 406 double y34 = (y3 + y4) / 2; 407 double x123 = (x12 + x23) / 2; 408 double y123 = (y12 + y23) / 2; 409 double x234 = (x23 + x34) / 2; 410 double y234 = (y23 + y34) / 2; 411 double x1234 = (x123 + x234) / 2; 412 double y1234 = (y123 + y234) / 2; 413 414 415 // Try to approximate the full cubic curve by a single straight line 416 //------------------ 417 double dx = x4-x1; 418 double dy = y4-y1; 419 420 double d2 = fabs(((x2 - x4) * dy - (y2 - y4) * dx)); 421 double d3 = fabs(((x3 - x4) * dy - (y3 - y4) * dx)); 422 double da1, da2, k; 423 424 switch((int(d2 > curve_collinearity_epsilon) << 1) + 425 int(d3 > curve_collinearity_epsilon)) 426 { 427 case 0: 428 // All collinear OR p1==p4 429 //---------------------- 430 k = dx*dx + dy*dy; 431 if(k == 0) 432 { 433 d2 = calc_sq_distance(x1, y1, x2, y2); 434 d3 = calc_sq_distance(x4, y4, x3, y3); 435 } 436 else 437 { 438 k = 1 / k; 439 da1 = x2 - x1; 440 da2 = y2 - y1; 441 d2 = k * (da1*dx + da2*dy); 442 da1 = x3 - x1; 443 da2 = y3 - y1; 444 d3 = k * (da1*dx + da2*dy); 445 if(d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1) 446 { 447 // Simple collinear case, 1---2---3---4 448 // We can leave just two endpoints 449 return; 450 } 451 if(d2 <= 0) d2 = calc_sq_distance(x2, y2, x1, y1); 452 else if(d2 >= 1) d2 = calc_sq_distance(x2, y2, x4, y4); 453 else d2 = calc_sq_distance(x2, y2, x1 + d2*dx, y1 + d2*dy); 454 455 if(d3 <= 0) d3 = calc_sq_distance(x3, y3, x1, y1); 456 else if(d3 >= 1) d3 = calc_sq_distance(x3, y3, x4, y4); 457 else d3 = calc_sq_distance(x3, y3, x1 + d3*dx, y1 + d3*dy); 458 } 459 if(d2 > d3) 460 { 461 if(d2 < m_distance_tolerance_square) 462 { 463 m_points.add(point_d(x2, y2)); 464 return; 465 } 466 } 467 else 468 { 469 if(d3 < m_distance_tolerance_square) 470 { 471 m_points.add(point_d(x3, y3)); 472 return; 473 } 474 } 475 break; 476 477 case 1: 478 // p1,p2,p4 are collinear, p3 is significant 479 //---------------------- 480 if(d3 * d3 <= m_distance_tolerance_square * (dx*dx + dy*dy)) 481 { 482 if(m_angle_tolerance < curve_angle_tolerance_epsilon) 483 { 484 m_points.add(point_d(x23, y23)); 485 return; 486 } 487 488 // Angle Condition 489 //---------------------- 490 da1 = fabs(atan2(y4 - y3, x4 - x3) - atan2(y3 - y2, x3 - x2)); 491 if(da1 >= pi) da1 = 2*pi - da1; 492 493 if(da1 < m_angle_tolerance) 494 { 495 m_points.add(point_d(x2, y2)); 496 m_points.add(point_d(x3, y3)); 497 return; 498 } 499 500 if(m_cusp_limit != 0.0) 501 { 502 if(da1 > m_cusp_limit) 503 { 504 m_points.add(point_d(x3, y3)); 505 return; 506 } 507 } 508 } 509 break; 510 511 case 2: 512 // p1,p3,p4 are collinear, p2 is significant 513 //---------------------- 514 if(d2 * d2 <= m_distance_tolerance_square * (dx*dx + dy*dy)) 515 { 516 if(m_angle_tolerance < curve_angle_tolerance_epsilon) 517 { 518 m_points.add(point_d(x23, y23)); 519 return; 520 } 521 522 // Angle Condition 523 //---------------------- 524 da1 = fabs(atan2(y3 - y2, x3 - x2) - atan2(y2 - y1, x2 - x1)); 525 if(da1 >= pi) da1 = 2*pi - da1; 526 527 if(da1 < m_angle_tolerance) 528 { 529 m_points.add(point_d(x2, y2)); 530 m_points.add(point_d(x3, y3)); 531 return; 532 } 533 534 if(m_cusp_limit != 0.0) 535 { 536 if(da1 > m_cusp_limit) 537 { 538 m_points.add(point_d(x2, y2)); 539 return; 540 } 541 } 542 } 543 break; 544 545 case 3: 546 // Regular case 547 //----------------- 548 if((d2 + d3)*(d2 + d3) <= m_distance_tolerance_square * (dx*dx + dy*dy)) 549 { 550 // If the curvature doesn't exceed the distance_tolerance value 551 // we tend to finish subdivisions. 552 //---------------------- 553 if(m_angle_tolerance < curve_angle_tolerance_epsilon) 554 { 555 m_points.add(point_d(x23, y23)); 556 return; 557 } 558 559 // Angle & Cusp Condition 560 //---------------------- 561 k = atan2(y3 - y2, x3 - x2); 562 da1 = fabs(k - atan2(y2 - y1, x2 - x1)); 563 da2 = fabs(atan2(y4 - y3, x4 - x3) - k); 564 if(da1 >= pi) da1 = 2*pi - da1; 565 if(da2 >= pi) da2 = 2*pi - da2; 566 567 if(da1 + da2 < m_angle_tolerance) 568 { 569 // Finally we can stop the recursion 570 //---------------------- 571 m_points.add(point_d(x23, y23)); 572 return; 573 } 574 575 if(m_cusp_limit != 0.0) 576 { 577 if(da1 > m_cusp_limit) 578 { 579 m_points.add(point_d(x2, y2)); 580 return; 581 } 582 583 if(da2 > m_cusp_limit) 584 { 585 m_points.add(point_d(x3, y3)); 586 return; 587 } 588 } 589 } 590 break; 591 } 592 593 // Continue subdivision 594 //---------------------- 595 recursive_bezier(x1, y1, x12, y12, x123, y123, x1234, y1234, level + 1); 596 recursive_bezier(x1234, y1234, x234, y234, x34, y34, x4, y4, level + 1); 597 } 598 599 //------------------------------------------------------------------------ 600 void curve4_div::bezier(double x1, double y1, 601 double x2, double y2, 602 double x3, double y3, 603 double x4, double y4) 604 { 605 m_points.add(point_d(x1, y1)); 606 recursive_bezier(x1, y1, x2, y2, x3, y3, x4, y4, 0); 607 m_points.add(point_d(x4, y4)); 608 } 609 610} 611 612