| Report problems to ATLAS LXR Team (with time and IP address indicated) |
|
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ] |
||||
|
||||||
| Links to LXR source navigation pages for stable releases | [ 12.*.* ] [ 13.*.* ] [ 14.*.* ] [ 15.*.* ] | |||||
001 #include "PolygonTriangulator.h" 002 003 /////////////////////////////////////////////////////////////////////// 004 // // 005 // Class: PolygonTriangulator // 006 // // 007 // Description: Triangulates any (non-complex) 2D polygon. // 008 // // 009 // Author: Thomas Kittelmann (Thomas.Kittelmann@cern.ch), March 2007 // 010 // // 011 // Notes: // 012 // // 013 // Actual algorithms are adapted from // 014 // http://www.mema.ucl.ac.be/~wu/Poly2Tri/poly2tri.html // 015 // (see copyright notice below) // 016 // // 017 // Main changes performed for the present incarnation are // 018 // interface & namespace changes, indentation, small // 019 // performance tweaks and removal of unused features. // 020 // Most importantly, got rid of static and globals so the class // 021 // can be reliably used more than once per process. Also, the // 022 // output triangles are sorted to preserve "handedness". // 023 // // 024 /////////////////////////////////////////////////////////////////////// 025 026 /////////////////////////////////////////////////////////////////////////////// 027 /////////////////////////////////////////////////////////////////////////////// 028 /////////////// Copyright file from poly2tri webpage follows: ///////////////// 029 /////////////////////////////////////////////////////////////////////////////// 030 /////////////////////////////////////////////////////////////////////////////// 031 // 032 // Poly2Tri:Fast and Robust Simple Polygon triangulation with/without holes 033 // by Sweep Line Algorithm 034 // Liang, Wu 035 // http://www.mema.ucl.ac.be/~wu/Poly2Tri/poly2tri.html 036 // Copyright (C) 2003, 2004, 2005, ALL RIGHTS RESERVED. 037 // 038 // --------------------------------------------------------------------- 039 // wu@mema.ucl.ac.be wuliang@femagsoft.com 040 // Centre for Sys. Eng. & App. Mech. FEMAGSoft S.A. 041 // Universite Cathalique de Louvain 4, Avenue Albert Einstein 042 // Batiment Euler, Avenue Georges Lemaitre, 4 B-1348 Louvain-la-Neuve 043 // B-1348, Louvain-la-Neuve Belgium 044 // Belgium 045 // --------------------------------------------------------------------- 046 // 047 // This program is distributed in the hope that it will be useful, but 048 // WITHOUT ANY WARRANTY; without even the implied warranty of MERCHAN- 049 // TABILITY or FITNESS FOR A PARTICULAR PURPOSE. 050 // 051 // This program may be freely redistributed under the condition that all 052 // the copyright notices in all source files ( including the copyright 053 // notice printed when the `-h' switch is selected) are not removed.Both 054 // the binary and source codes may not be sold or included in any comme- 055 // rcial products without a license from the corresponding author(s) & 056 // entities. 057 // 058 // 1) Arbitrary precision floating-point arithmetic and fast robust geo- 059 // metric predicates (predicates.cc) is copyrighted by 060 // Jonathan Shewchuk (http://www.cs.berkeley.edu/~jrs) and you may get 061 // the source code from http://www.cs.cmu.edu/~quake/robust.html 062 // 063 // 2) The shell script mps2eps is copyrighted by Jon Edvardsson 064 // (http://www.ida.liu.se/~pelab/members/index.php4/?12) and you may 065 // get the copy from http://www.ida.liu.se/~joned/download/mps2eps/ 066 // 067 // 3) All other source codes and exmaples files distributed in Poly2Tri 068 // are copyrighted by Liang, Wu (http://www.mema.ucl.ac.be/~wu) and 069 // FEMAGSoft S.A. 070 // 071 /////////////////////////////////////////////////////////////////////////////// 072 /////////////////////////////////////////////////////////////////////////////// 073 /////////////////////////////////////////////////////////////////////////////// 074 075 #include <algorithm> 076 #include <cmath> 077 #include <cstdlib> 078 #include <map> 079 #include <limits> 080 #include <list> 081 #include <queue> 082 #include <set> 083 #include <iostream> 084 #include <stack> 085 #include <assert.h> 086 087 namespace internal_poltrig { 088 089 //////////////////////////////////////////////////////////////////////////////////////////// 090 //////////////////////////////////////////////////////////////////////////////////////////// 091 //////////////////////////////////// Defs.h //////////////////////////////////////////////// 092 //////////////////////////////////////////////////////////////////////////////////////////// 093 //////////////////////////////////////////////////////////////////////////////////////////// 094 095 //-------------------------------------------------------------------------/ 096 //Copyright (C) 2003, 2004, 2005, ALL RIGHTS RESERVED. 097 //Centre for Sys. Eng. & App. Mech. FEMAGSoft S.A. 098 //Universite Cathalique de Louvain 4, Avenue Albert Einstein 099 //Batiment Euler, Avenue Georges Lemaitre, 4 B-1348 Louvain-la-Neuve 100 //B-1348, Louvain-la-Neuve Belgium 101 //Belgium 102 //-------------------------------------------------------------------------/ 103 //Name: defs.h 104 //Purpose: some definition for mesh generation 105 //Author: Wu Liang 106 //Created: 03/2000 107 //-------------------------------------------------------------------------/ 108 109 110 #define sqr(t) (t)*(t) 111 112 const double PI=3.141592653589793238462643383279502884197169399375105820974944592308; 113 enum Type { UNKNOWN, INPUT, INSERT, START, END, MERGE, SPLIT, REGULAR_UP, REGULAR_DOWN}; 114 115 class Pointbase; 116 class Linebase; 117 118 template <class T, class KeyType> class SplayTree; 119 typedef std::map<unsigned int, Pointbase*> PointbaseMap; 120 typedef std::map<unsigned int, Linebase*> LineMap; 121 typedef std::priority_queue<Pointbase> PQueue; 122 typedef SplayTree<Linebase*, double> EdgeBST; 123 typedef std::list<unsigned int> Monopoly; 124 typedef std::list<Monopoly> Monopolys; 125 typedef std::vector<unsigned int> Triangle; 126 typedef std::list<Triangle> Triangles; 127 typedef std::map<unsigned int, std::set<unsigned int> > AdjEdgeMap; 128 129 130 //////////////////////////////////////////////////////////////////////////////////////////// 131 //////////////////////////////////////////////////////////////////////////////////////////// 132 //////////////////////////////// predicates.cc ///////////////////////////////////////////// 133 //////////////////////////////////////////////////////////////////////////////////////////// 134 //////////////////////////////////////////////////////////////////////////////////////////// 135 136 /*****************************************************************************/ 137 /* */ 138 /* Routines for Arbitrary Precision Floating-point Arithmetic */ 139 /* and Fast Robust Geometric Predicates */ 140 /* (predicates.cc) */ 141 /* */ 142 /* May 18, 1996 */ 143 /* */ 144 /* Placed in the public domain by */ 145 /* Jonathan Richard Shewchuk */ 146 /* School of Computer Science */ 147 /* Carnegie Mellon University */ 148 /* 5000 Forbes Avenue */ 149 /* Pittsburgh, Pennsylvania 15213-3891 */ 150 /* jrs@cs.cmu.edu */ 151 /* */ 152 /* This file contains C implementation of algorithms for exact addition */ 153 /* and multiplication of floating-point numbers, and predicates for */ 154 /* robustly performing the orientation and incircle tests used in */ 155 /* computational geometry. The algorithms and underlying theory are */ 156 /* described in Jonathan Richard Shewchuk. "Adaptive Precision Floating- */ 157 /* Point Arithmetic and Fast Robust Geometric Predicates." Technical */ 158 /* Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon */ 159 /* University, Pittsburgh, Pennsylvania, May 1996. (Submitted to */ 160 /* Discrete & Computational Geometry.) */ 161 /* */ 162 /* This file, the paper listed above, and other information are available */ 163 /* from the Web page http://www.cs.cmu.edu/~quake/robust.html . */ 164 /* */ 165 /*****************************************************************************/ 166 167 /*****************************************************************************/ 168 /* */ 169 /* Using this code: */ 170 /* */ 171 /* First, read the short or long version of the paper (from the Web page */ 172 /* above). */ 173 /* */ 174 /* Be sure to call exactinit() once, before calling any of the arithmetic */ 175 /* functions or geometric predicates. Also be sure to turn on the */ 176 /* optimizer when compiling this file. */ 177 /* */ 178 /* */ 179 /* Several geometric predicates are defined. Their parameters are all */ 180 /* points. Each point is an array of two or three floating-point */ 181 /* numbers. The geometric predicates, described in the papers, are */ 182 /* */ 183 /* orient2d(pa, pb, pc) */ 184 /* orient2dfast(pa, pb, pc) */ 185 /* orient3d(pa, pb, pc, pd) */ 186 /* orient3dfast(pa, pb, pc, pd) */ 187 /* incircle(pa, pb, pc, pd) */ 188 /* incirclefast(pa, pb, pc, pd) */ 189 /* insphere(pa, pb, pc, pd, pe) */ 190 /* inspherefast(pa, pb, pc, pd, pe) */ 191 /* */ 192 /* Those with suffix "fast" are approximate, non-robust versions. Those */ 193 /* without the suffix are adaptive precision, robust versions. There */ 194 /* are also versions with the suffices "exact" and "slow", which are */ 195 /* non-adaptive, exact arithmetic versions, which I use only for timings */ 196 /* in my arithmetic papers. */ 197 /* */ 198 /* */ 199 /* An expansion is represented by an array of floating-point numbers, */ 200 /* sorted from smallest to largest magnitude (possibly with interspersed */ 201 /* zeros). The length of each expansion is stored as a separate integer, */ 202 /* and each arithmetic function returns an integer which is the length */ 203 /* of the expansion it created. */ 204 /* */ 205 /* Several arithmetic functions are defined. Their parameters are */ 206 /* */ 207 /* e, f Input expansions */ 208 /* elen, flen Lengths of input expansions (must be >= 1) */ 209 /* h Output expansion */ 210 /* b Input scalar */ 211 /* */ 212 /* The arithmetic functions are */ 213 /* */ 214 /* grow_expansion(elen, e, b, h) */ 215 /* grow_expansion_zeroelim(elen, e, b, h) */ 216 /* expansion_sum(elen, e, flen, f, h) */ 217 /* expansion_sum_zeroelim1(elen, e, flen, f, h) */ 218 /* expansion_sum_zeroelim2(elen, e, flen, f, h) */ 219 /* fast_expansion_sum(elen, e, flen, f, h) */ 220 /* fast_expansion_sum_zeroelim(elen, e, flen, f, h) */ 221 /* linear_expansion_sum(elen, e, flen, f, h) */ 222 /* linear_expansion_sum_zeroelim(elen, e, flen, f, h) */ 223 /* scale_expansion(elen, e, b, h) */ 224 /* scale_expansion_zeroelim(elen, e, b, h) */ 225 /* compress(elen, e, h) */ 226 /* */ 227 /* All of these are described in the long version of the paper; some are */ 228 /* described in the short version. All return an integer that is the */ 229 /* length of h. Those with suffix _zeroelim perform zero elimination, */ 230 /* and are recommended over their counterparts. The procedure */ 231 /* fast_expansion_sum_zeroelim() (or linear_expansion_sum_zeroelim() on */ 232 /* processors that do not use the round-to-even tiebreaking rule) is */ 233 /* recommended over expansion_sum_zeroelim(). Each procedure has a */ 234 /* little note next to it (in the code below) that tells you whether or */ 235 /* not the output expansion may be the same array as one of the input */ 236 /* expansions. */ 237 /* */ 238 /* */ 239 /* If you look around below, you'll also find macros for a bunch of */ 240 /* simple unrolled arithmetic operations, and procedures for printing */ 241 /* expansions (commented out because they don't work with all C */ 242 /* compilers) and for generating rand floating-point numbers whose */ 243 /* significand bits are all rand. Most of the macros have undocumented */ 244 /* requirements that certain of their parameters should not be the same */ 245 /* variable; for safety, better to make sure all the parameters are */ 246 /* distinct variables. Feel free to send email to jrs@cs.cmu.edu if you */ 247 /* have questions. */ 248 /* */ 249 /*****************************************************************************/ 250 251 /* On some machines, the exact arithmetic routines might be defeated by the */ 252 /* use of internal extended precision floating-point registers. Sometimes */ 253 /* this problem can be fixed by defining certain values to be volatile, */ 254 /* thus forcing them to be stored to memory and rounded off. This isn't */ 255 /* a great solution, though, as it slows the arithmetic down. */ 256 /* */ 257 /* To try this out, write "#define INEXACT volatile" below. Normally, */ 258 /* however, INEXACT should be defined to be nothing. ("#define INEXACT".) */ 259 260 #define INEXACT /* Nothing */ 261 /* #define INEXACT volatile */ 262 263 #define REAL double /* float or double */ 264 265 /* Which of the following two methods of finding the absolute values is */ 266 /* fastest is compiler-dependent. A few compilers can inline and optimize */ 267 /* the fabs() call; but most will incur the overhead of a function call, */ 268 /* which is disastrously slow. A faster way on IEEE machines might be to */ 269 /* mask the appropriate bit, but that's difficult to do in C. */ 270 271 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a)) 272 /* #define Absolute(a) fabs(a) */ 273 274 /* Many of the operations are broken up into two pieces, a main part that */ 275 /* performs an approximate operation, and a "tail" that computes the */ 276 /* roundoff error of that operation. */ 277 /* */ 278 /* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), */ 279 /* Split(), and Two_Product() are all implemented as described in the */ 280 /* reference. Each of these macros requires certain variables to be */ 281 /* defined in the calling routine. The variables `bvirt', `c', `abig', */ 282 /* `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because */ 283 /* they store the result of an operation that may incur roundoff error. */ 284 /* The input parameter `x' (or the highest numbered `x_' parameter) must */ 285 /* also be declared `INEXACT'. */ 286 287 #define Fast_Two_Sum_Tail(a, b, x, y) \ 288 bvirt = x - a; \ 289 y = b - bvirt 290 291 #define Fast_Two_Sum(a, b, x, y) \ 292 x = (REAL) (a + b); \ 293 Fast_Two_Sum_Tail(a, b, x, y) 294 295 #define Two_Sum_Tail(a, b, x, y) \ 296 bvirt = (REAL) (x - a); \ 297 avirt = x - bvirt; \ 298 bround = b - bvirt; \ 299 around = a - avirt; \ 300 y = around + bround 301 302 #define Two_Sum(a, b, x, y) \ 303 x = (REAL) (a + b); \ 304 Two_Sum_Tail(a, b, x, y) 305 306 #define Two_Diff_Tail(a, b, x, y) \ 307 bvirt = (REAL) (a - x); \ 308 avirt = x + bvirt; \ 309 bround = bvirt - b; \ 310 around = a - avirt; \ 311 y = around + bround 312 313 #define Two_Diff(a, b, x, y) \ 314 x = (REAL) (a - b); \ 315 Two_Diff_Tail(a, b, x, y) 316 317 318 // c = (REAL) (splitter * a); 319 320 #define Split(a, ahi, alo) \ 321 c = (REAL) (1.0 * a); \ 322 abig = (REAL) (c - a); \ 323 ahi = c - abig; \ 324 alo = a - ahi 325 326 #define Two_Product_Tail(a, b, x, y) \ 327 Split(a, ahi, alo); \ 328 Split(b, bhi, blo); \ 329 err1 = x - (ahi * bhi); \ 330 err2 = err1 - (alo * bhi); \ 331 err3 = err2 - (ahi * blo); \ 332 y = (alo * blo) - err3 333 334 #define Two_Product(a, b, x, y) \ 335 x = (REAL) (a * b); \ 336 Two_Product_Tail(a, b, x, y) 337 338 /* Macros for summing expansions of various fixed lengths. These are all */ 339 /* unrolled versions of Expansion_Sum(). */ 340 341 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \ 342 Two_Diff(a0, b , _i, x0); \ 343 Two_Sum( a1, _i, x2, x1) 344 345 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \ 346 Two_One_Diff(a1, a0, b0, _j, _0, x0); \ 347 Two_One_Diff(_j, _0, b1, x3, x2, x1) 348 349 350 //TK: weird globals, never initialised. Just replaced with 1.0 everywhere... 351 // REAL splitter(1); /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */ 352 /* A set of coefficients used to calculate maximum roundoff errors. */ 353 // REAL resulterrbound; 354 // REAL ccwerrboundA, ccwerrboundB, ccwerrboundC; 355 // REAL o3derrboundA, o3derrboundB, o3derrboundC; 356 // REAL iccerrboundA, iccerrboundB, iccerrboundC; 357 // REAL isperrboundA, isperrboundB, isperrboundC; 358 359 /*****************************************************************************/ 360 /* */ 361 /* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */ 362 /* components from the output expansion. */ 363 /* */ 364 /* Sets h = e + f. See the long version of my paper for details. */ 365 /* */ 366 /* If round-to-even is used (as with IEEE 754), maintains the strongly */ 367 /* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */ 368 /* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */ 369 /* properties. */ 370 /* */ 371 /*****************************************************************************/ 372 373 int fast_expansion_sum_zeroelim(const int& elen, REAL* e, const int& flen, REAL* f, REAL* h) 374 /* h cannot be e or f. */ 375 { 376 REAL Q; 377 INEXACT REAL Qnew; 378 INEXACT REAL hh; 379 INEXACT REAL bvirt; 380 REAL avirt, bround, around; 381 int eindex, findex, hindex; 382 REAL enow, fnow; 383 384 enow = e[0]; 385 fnow = f[0]; 386 eindex = findex = 0; 387 if ((fnow > enow) == (fnow > -enow)) { 388 Q = enow; 389 enow = e[++eindex]; 390 } else { 391 Q = fnow; 392 fnow = f[++findex]; 393 } 394 hindex = 0; 395 if ((eindex < elen) && (findex < flen)) { 396 if ((fnow > enow) == (fnow > -enow)) { 397 Fast_Two_Sum(enow, Q, Qnew, hh); 398 enow = e[++eindex]; 399 } else { 400 Fast_Two_Sum(fnow, Q, Qnew, hh); 401 fnow = f[++findex]; 402 } 403 Q = Qnew; 404 if (hh != 0.0) { 405 h[hindex++] = hh; 406 } 407 while ((eindex < elen) && (findex < flen)) { 408 if ((fnow > enow) == (fnow > -enow)) { 409 Two_Sum(Q, enow, Qnew, hh); 410 enow = e[++eindex]; 411 } else { 412 Two_Sum(Q, fnow, Qnew, hh); 413 fnow = f[++findex]; 414 } 415 Q = Qnew; 416 if (hh != 0.0) { 417 h[hindex++] = hh; 418 } 419 } 420 } 421 while (eindex < elen) { 422 Two_Sum(Q, enow, Qnew, hh); 423 enow = e[++eindex]; 424 Q = Qnew; 425 if (hh != 0.0) { 426 h[hindex++] = hh; 427 } 428 } 429 while (findex < flen) { 430 Two_Sum(Q, fnow, Qnew, hh); 431 fnow = f[++findex]; 432 Q = Qnew; 433 if (hh != 0.0) { 434 h[hindex++] = hh; 435 } 436 } 437 if ((Q != 0.0) || (hindex == 0)) { 438 h[hindex++] = Q; 439 } 440 return hindex; 441 } 442 443 444 /*****************************************************************************/ 445 /* */ 446 /* estimate() Produce a one-word estimate of an expansion's value. */ 447 /* */ 448 /* See either version of my paper for details. */ 449 /* */ 450 /*****************************************************************************/ 451 452 REAL estimate(const int& elen, REAL* e) 453 { 454 REAL Q; 455 int eindex; 456 457 Q = e[0]; 458 for (eindex = 1; eindex < elen; ++eindex) { 459 Q += e[eindex]; 460 } 461 return Q; 462 } 463 464 /*****************************************************************************/ 465 /* */ 466 /* orient2dfast() Approximate 2D orientation test. Nonrobust. */ 467 /* orient2dexact() Exact 2D orientation test. Robust. */ 468 /* orient2dslow() Another exact 2D orientation test. Robust. */ 469 /* orient2d() Adaptive exact 2D orientation test. Robust. */ 470 /* */ 471 /* Return a positive value if the points pa, pb, and pc occur */ 472 /* in counterclockwise order; a negative value if they occur */ 473 /* in clockwise order; and zero if they are collinear. The */ 474 /* result is also a rough approximation of twice the signed */ 475 /* area of the triangle defined by the three points. */ 476 /* */ 477 /* Only the first and last routine should be used; the middle two are for */ 478 /* timings. */ 479 /* */ 480 /* The last three use exact arithmetic to ensure a correct answer. The */ 481 /* result returned is the determinant of a matrix. In orient2d() only, */ 482 /* this determinant is computed adaptively, in the sense that exact */ 483 /* arithmetic is used only to the degree it is needed to ensure that the */ 484 /* returned value has the correct sign. Hence, orient2d() is usually quite */ 485 /* fast, but will run more slowly when the input points are collinear or */ 486 /* nearly so. */ 487 /* */ 488 /*****************************************************************************/ 489 490 REAL orient2dadapt(REAL* pa, REAL* pb, REAL* pc, const REAL& detsum) 491 { 492 INEXACT REAL acx, acy, bcx, bcy; 493 REAL acxtail, acytail, bcxtail, bcytail; 494 INEXACT REAL detleft, detright; 495 REAL detlefttail, detrighttail; 496 REAL det, errbound; 497 REAL B[4], C1[8], C2[12], D[16]; 498 INEXACT REAL B3; 499 int C1length, C2length, Dlength; 500 REAL u[4]; 501 INEXACT REAL u3; 502 INEXACT REAL s1, t1; 503 REAL s0, t0; 504 505 INEXACT REAL bvirt; 506 REAL avirt, bround, around; 507 INEXACT REAL c; 508 INEXACT REAL abig; 509 REAL ahi, alo, bhi, blo; 510 REAL err1, err2, err3; 511 INEXACT REAL _i, _j; 512 REAL _0; 513 514 acx = (REAL) (pa[0] - pc[0]); 515 bcx = (REAL) (pb[0] - pc[0]); 516 acy = (REAL) (pa[1] - pc[1]); 517 bcy = (REAL) (pb[1] - pc[1]); 518 519 Two_Product(acx, bcy, detleft, detlefttail); 520 Two_Product(acy, bcx, detright, detrighttail); 521 522 Two_Two_Diff(detleft, detlefttail, detright, detrighttail, 523 B3, B[2], B[1], B[0]); 524 B[3] = B3; 525 526 det = estimate(4, B); 527 // errbound = ccwerrboundB * detsum; 528 errbound = detsum; 529 if ((det >= errbound) || (-det >= errbound)) { 530 return det; 531 } 532 533 Two_Diff_Tail(pa[0], pc[0], acx, acxtail); 534 Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail); 535 Two_Diff_Tail(pa[1], pc[1], acy, acytail); 536 Two_Diff_Tail(pb[1], pc[1], bcy, bcytail); 537 538 if ((acxtail == 0.0) && (acytail == 0.0) 539 && (bcxtail == 0.0) && (bcytail == 0.0)) { 540 return det; 541 } 542 543 // errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det); 544 errbound = detsum + Absolute(det); 545 det += (acx * bcytail + bcy * acxtail) 546 - (acy * bcxtail + bcx * acytail); 547 if ((det >= errbound) || (-det >= errbound)) { 548 return det; 549 } 550 551 Two_Product(acxtail, bcy, s1, s0); 552 Two_Product(acytail, bcx, t1, t0); 553 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); 554 u[3] = u3; 555 C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1); 556 557 Two_Product(acx, bcytail, s1, s0); 558 Two_Product(acy, bcxtail, t1, t0); 559 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); 560 u[3] = u3; 561 C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2); 562 563 Two_Product(acxtail, bcytail, s1, s0); 564 Two_Product(acytail, bcxtail, t1, t0); 565 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); 566 u[3] = u3; 567 Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D); 568 569 return(D[Dlength - 1]); 570 } 571 572 REAL orient2d(REAL* pa, REAL* pb, REAL* pc) 573 { 574 REAL detleft, detright, det; 575 REAL detsum, errbound; 576 577 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]); 578 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]); 579 det = detleft - detright; 580 581 if (detleft > 0.0) { 582 if (detright <= 0.0) { 583 return det; 584 } else { 585 detsum = detleft + detright; 586 } 587 } else if (detleft < 0.0) { 588 if (detright >= 0.0) { 589 return det; 590 } else { 591 detsum = -detleft - detright; 592 } 593 } else { 594 return det; 595 } 596 597 // errbound = ccwerrboundA * detsum; 598 errbound = detsum; 599 if ((det >= errbound) || (-det >= errbound)) { 600 return det; 601 } 602 603 return orient2dadapt(pa, pb, pc, detsum); 604 } 605 606 607 //////////////////////////////////////////////// 608 609 //////////////////////////////////////////////////////////////////////////////////////////// 610 //////////////////////////////////////////////////////////////////////////////////////////// 611 //////////////////////////////////// splay.h /////////////////////////////////////////////// 612 //////////////////////////////////////////////////////////////////////////////////////////// 613 //////////////////////////////////////////////////////////////////////////////////////////// 614 615 //-------------------------------------------------------------------------/ 616 //Copyright (C) 2003, 2004, 2005, ALL RIGHTS RESERVED. 617 //Centre for Sys. Eng. & App. Mech. FEMAGSoft S.A. 618 //Universite Cathalique de Louvain 4, Avenue Albert Einstein 619 //Batiment Euler, Avenue Georges Lemaitre, 4 B-1348 Louvain-la-Neuve 620 //B-1348, Louvain-la-Neuve Belgium 621 //Belgium 622 //-------------------------------------------------------------------------/ 623 // 624 //Name: splay.h (self-adjusting balanced binary search tree) 625 //Purpose: declaration and implentation of splay tree class for fast 626 // binary searching 627 //Author: Liang, Wu (wu@mema.ucl.ac.be, wuliang@femagsoft.com) 628 //Created: 03/2002 629 //Modified: 2003, 2004, 2005 630 //-------------------------------------------------------------------------/ 631 632 template <class T, class KeyType> 633 class SplayTree; 634 635 template <class T, class KeyType> 636 class BTreeNode 637 { 638 public: 639 friend class SplayTree<T, KeyType>; 640 BTreeNode( ) : _left( NULL ), _right( NULL ), _visited(false) { } 641 BTreeNode( const T & data, BTreeNode *lt, BTreeNode *rt ) 642 : _data(data),_left( lt ), _right( rt ), _visited(false) { } 643 644 T& data() { return _data; } 645 BTreeNode* Left() { return _left; } 646 BTreeNode* Right() { return _right; } 647 void SetVisited(const bool& visited) { _visited=visited; } 648 bool GetVisited() { return _visited; } 649 KeyType keyValue() { return _data->keyValue(); } 650 651 private: 652 T _data; 653 BTreeNode *_left; 654 BTreeNode *_right; 655 bool _visited; 656 657 }; 658 659 660 template <class T, class KeyType> 661 class SplayTree 662 { 663 public: 664 explicit SplayTree( ):root(NULL),size(0) { } 665 SplayTree( const SplayTree & rhs ); 666 ~SplayTree( ); 667 668 void MakeEmpty( ); 669 bool IsEmpty( ) const; 670 long int Size() { return size; } 671 BTreeNode<T, KeyType>* Root() { return root; } 672 673 void Find( const KeyType& keys, BTreeNode<T, KeyType>* & res); 674 void FindMin( BTreeNode<T, KeyType>* &min ); 675 void FindMax( BTreeNode<T, KeyType>* &max ); 676 // We only use this function for polygon Triangulation to find the direct 677 // left edge of an event vertex; 678 void FindMaxSmallerThan( const KeyType& keys, BTreeNode<T, KeyType>* &res); 679 680 void Insert( const T & x ); 681 void Delete( const KeyType& keys); 682 void Delete( const KeyType& keys, BTreeNode<T, KeyType>* &res); 683 void DeleteMin(BTreeNode<T, KeyType>* &min); 684 void DeleteMax(BTreeNode<T, KeyType>* &max); 685 686 const SplayTree & operator=( const SplayTree & rhs ); 687 void PreOrder( void(*Visit)(BTreeNode<T,KeyType> *u) ) 688 { PreOrder(Visit, root); } 689 void InOrder( void(*Visit)(BTreeNode<T,KeyType> *u) ) 690 { InOrder(Visit, root); } 691 692 void InOrder( void(*Visit)(BTreeNode<T,KeyType>*u, double y), double y) 693 { InOrder(Visit, root, y); } 694 695 void PostOrder( void(*Visit)(BTreeNode<T,KeyType> *u) ) 696 { PostOrder(Visit, root); } 697 698 int Height( ) const { return Height(root); } //height of root 699 int Height(BTreeNode<T, KeyType> *t) const; //Height of subtree t; 700 BTreeNode<T, KeyType>* Left(BTreeNode<T, KeyType> *node) { return node->_left; } 701 BTreeNode<T, KeyType>* Right(BTreeNode<T, KeyType> *node) { return node->_right; } 702 703 private: 704 BTreeNode<T, KeyType> *root; 705 long int size; 706 707 void reclaimMemory( BTreeNode<T, KeyType> * t ) const; 708 BTreeNode<T, KeyType> * clone( BTreeNode<T, KeyType> *t ) const; 709 710 //Transverse 711 void PreOrder( void(*Visit)(BTreeNode<T, KeyType> *u), BTreeNode<T, KeyType> *t); 712 void InOrder( void(*Visit)(BTreeNode<T, KeyType> *u), BTreeNode<T, KeyType> *t); 713 void PostOrder( void(*Visit)(BTreeNode<T, KeyType> *u), BTreeNode<T, KeyType> *t); 714 715 void InOrder( void(*Visit)(BTreeNode<T, KeyType>*, double y), 716 BTreeNode<T, KeyType> *t, double y); 717 718 719 720 // Tree manipulations 721 void rotateWithLeftChild( BTreeNode<T, KeyType> * & k2 ) const; 722 void rotateWithRightChild( BTreeNode<T, KeyType> * & k1 ) const; 723 void splay( const KeyType& keys, BTreeNode<T, KeyType> * & t ) const; 724 }; 725 726 //---------------------------------------------------------------------- 727 //Constructor; 728 //---------------------------------------------------------------------- 729 template <class T, class KeyType> 730 SplayTree<T, KeyType>::SplayTree( const SplayTree<T, KeyType> & rhs ) 731 { 732 *this = rhs; 733 } 734 735 //----------------------------------------------------------------------- 736 //Destructor. 737 //----------------------------------------------------------------------- 738 template <class T, class KeyType> 739 SplayTree<T, KeyType>::~SplayTree( ) 740 { 741 MakeEmpty( ); 742 } 743 744 //------------------------------------------------------------------------ 745 //Insert x into the tree. 746 //------------------------------------------------------------------------ 747 template <class T, class KeyType> 748 void SplayTree<T, KeyType>::Insert( const T & x ) 749 { 750 751 BTreeNode<T, KeyType> *newNode= new BTreeNode<T, KeyType>; 752 newNode->_data=x; 753 754 if( root == NULL ) 755 { 756 newNode->_left = newNode->_right = NULL; 757 root = newNode; ++size; 758 } 759 else 760 { 761 KeyType keys=x->keyValue(); 762 splay( keys, root ); 763 KeyType rootk=root->keyValue(); 764 if( keys < rootk ) 765 { 766 newNode->_left = root->_left; 767 newNode->_right = root; 768 root->_left = NULL; 769 root = newNode; 770 ++size; 771 } 772 else if( keys > rootk ) 773 { 774 775 newNode->_right = root->_right; 776 newNode->_left = root; 777 root->_right = NULL; 778 root = newNode; 779 ++size; 780 } 781 else 782 { 783 //slight incresed the keyvalue to avoid duplicated keys 784 x->increaseKeyValue(1.0e-10); 785 Insert(x); 786 787 } 788 } 789 } 790 791 //--------------------------------------------------------------------- 792 //Remove the node with the keys from the tree, and retrieve the result 793 //--------------------------------------------------------------------- 794 template <class T, class KeyType> 795 void SplayTree<T, KeyType>::Delete( const KeyType& keys, BTreeNode<T, KeyType>* &res) 796 { 797 BTreeNode<T, KeyType> *newTree; 798 799 splay( keys, root ); 800 if( root->keyValue() != keys ) { res=NULL; return; } // Item not found; do nothing 801 802 res = root; 803 804 if( root->_left == NULL ) 805 newTree = root->_right; 806 else 807 { 808 // Find the maximum in the _left subtree 809 // Splay it to the root; and then attach _right child 810 newTree = root->_left; 811 splay( keys, newTree ); 812 newTree->_right = root->_right; 813 } 814 815 root = newTree; 816 size--; 817 } 818 819 //--------------------------------------------------------------------- 820 //Remove the node with the keys from the tree. 821 //--------------------------------------------------------------------- 822 template <class T, class KeyType> 823 void SplayTree<T, KeyType>::Delete( const KeyType& keys) 824 { 825 BTreeNode<T, KeyType> *newTree; 826 827 splay( keys, root ); 828 KeyType rootk=root->keyValue(); 829 if( rootk != keys ) { return; } // Item not found; do nothing 830 831 if( root->_left == NULL ) newTree = root->_right; 832 else 833 { 834 // Find the maximum in the _left subtree 835 // Splay it to the root; and then attach _right child 836 newTree = root->_left; 837 splay( keys, newTree ); 838 newTree->_right = root->_right; 839 } 840 841 delete root; 842 root = newTree; 843 size--; 844 } 845 846 847 848 //--------------------------------------------------------------------- 849 //Find and Delete the node with min keys from the tree; 850 //--------------------------------------------------------------------- 851 template <class T, class KeyType> 852 void SplayTree<T, KeyType>::DeleteMin(BTreeNode<T, KeyType>* &min) 853 { 854 if( IsEmpty( ) ) { min=NULL; return; } 855 856 double keys=-1.0e30; 857 splay( keys, root ); 858 859 min = root; 860 861 BTreeNode<T, KeyType> *newTree; 862 if( root->_left == NULL ) newTree = root->_right; 863 else 864 { 865 newTree = root->_left; 866 splay( keys, newTree ); 867 newTree->_right = root->_right; 868 } 869 870 size--; 871 root = newTree; 872 873 } 874 875 //---------------------------------------------------------------------- 876 //Find and Delete the node with max keys from the tree; 877 //---------------------------------------------------------------------- 878 template <class T, class KeyType> 879 void SplayTree<T, KeyType>::DeleteMax(BTreeNode<T, KeyType>* &max) 880 { 881 if( IsEmpty( ) ) { max=NULL; return; } 882 883 double keys=1.0e30; 884 splay( keys, root ); 885 886 max = root; 887 888 BTreeNode<T, KeyType> *newTree; 889 if( root->_left == NULL ) newTree = root->_right; 890 else 891 { 892 newTree = root->_left; 893 splay( keys, newTree ); 894 newTree->_right = root->_right; 895 } 896 size--; 897 root = newTree; 898 } 899 900 901 //---------------------------------------------------------------------- 902 //Find the smallest item in the tree, but won't delete it from the tree. 903 //---------------------------------------------------------------------- 904 template <class T, class KeyType> 905 void SplayTree<T, KeyType>::FindMin(BTreeNode<T, KeyType>* & min ) 906 { 907 if( IsEmpty( ) ) { min=NULL; return; } 908 BTreeNode<T, KeyType> *ptr = root; 909 910 while( ptr->_left != NULL ) ptr = ptr->_left; 911 splay( ptr->keyValue(), root ); 912 min = ptr; 913 } 914 915 //---------------------------------------------------------------------- 916 //Find the largest item in the tree. but won't delete it from the tree. 917 //---------------------------------------------------------------------- 918 template <class T, class KeyType> 919 void SplayTree<T, KeyType>::FindMax(BTreeNode<T, KeyType>* & max) 920 { 921 if( IsEmpty( ) ) { max=NULL; return; } 922 923 BTreeNode<T, KeyType> *ptr = root; 924 while( ptr->_right != NULL ) ptr = ptr->_right; 925 splay( ptr->keyValue(), root ); 926 max = ptr; 927 } 928 929 //-------------------------------------------------------------------- 930 //Find the node with the keys in the tree. 931 //res==NULL if it donesn't exist in the tree; 932 //-------------------------------------------------------------------- 933 template <class T, class KeyType> 934 void SplayTree<T, KeyType>::Find( const KeyType& keys, BTreeNode<T, KeyType>* & res) 935 { 936 if( IsEmpty( ) ) { res=NULL; return; } 937 splay( keys, root ); 938 if( root->keyValue() != keys ) { res=NULL; return; } 939 else res = root; 940 } 941 942 //-------------------------------------------------------------------- 943 //Find the maximum node smaller than or equal to the given key. 944 //This function specially designed for polygon Triangulation to 945 //find the direct left edge at event vertex; 946 //-------------------------------------------------------------------- 947 template <class T, class KeyType> 948 void SplayTree<T, KeyType>::FindMaxSmallerThan( const KeyType& keys, BTreeNode<T, KeyType>* &res) 949 { 950 if( IsEmpty( ) ) { res=NULL; return; } 951 splay( keys, root ); 952 953 if( root->data()->keyValue() < keys) res=root; 954 else if(root->_left) 955 { 956 res=root->_left; 957 while(res->_right) res=res->_right; 958 } 959 else 960 { 961 assert(false); 962 } 963 } 964 965 //-------------------------------------------------------------------- 966 //Make the tree logically empty. 967 //-------------------------------------------------------------------- 968 template <class T, class KeyType> 969 void SplayTree<T, KeyType>::MakeEmpty( ) 970 { 971 BTreeNode<T, KeyType>* ptr; 972 while( !IsEmpty( ) ) 973 { 974 DeleteMax(ptr); 975 delete ptr; 976 } 977 } 978 979 //--------------------------------------------------------------------- 980 //Test if the tree is logically empty. 981 //--------------------------------------------------------------------- 982 template <class T, class KeyType> 983 bool SplayTree<T, KeyType>::IsEmpty( ) const 984 { 985 return root == NULL; 986 } 987 988 //---------------------------------------------------------------------- 989 //copy overload operator. 990 //---------------------------------------------------------------------- 991 template <class T, class KeyType> 992 const SplayTree<T, KeyType> & SplayTree<T, KeyType>::operator=( const SplayTree<T, KeyType> & rhs ) 993 { 994 if( this != &rhs ) 995 { 996 MakeEmpty( ); 997 root = clone( rhs.root ); 998 } 999 1000 return *this; 1001 } 1002 1003 //----------------------------------------------------------------------- 1004 //Internal method to perform a top-down splay. 1005 //x is the key of target node to splay around. 1006 //t is the root of the subtree to splay. 1007 //----------------------------------------------------------------------- 1008 template <class T, class KeyType> 1009 void SplayTree<T, KeyType>::splay( const KeyType& keys, BTreeNode<T, KeyType> * & t ) const 1010 { 1011 BTreeNode<T, KeyType> *_leftTreeMax, *_rightTreeMin; 1012 // static BTreeNode<T, KeyType> header; 1013 BTreeNode<T, KeyType> header;//TK: Removed static keyword. Rather a bit slower than thread problems... 1014 1015 header._left = header._right = NULL; 1016 _leftTreeMax = _rightTreeMin = &header; 1017 1018 for( ; ; ) 1019 { 1020 KeyType rkey=t->keyValue(); 1021 if( keys < rkey ) 1022 { 1023 if(t->_left == NULL) break; 1024 if( keys < t->_left->keyValue() ) rotateWithLeftChild( t ); 1025 if( t->_left == NULL ) break; 1026 1027 // Link Right 1028 _rightTreeMin->_left = t; 1029 _rightTreeMin = t; 1030 t = t->_left; 1031 } 1032 else if( keys > rkey ) 1033 { 1034 if( t->_right == NULL ) break; 1035 if( keys > t->_right->keyValue() ) rotateWithRightChild( t ); 1036 if( t->_right == NULL ) break; 1037 1038 // Link Left 1039 _leftTreeMax->_right = t; 1040 _leftTreeMax = t; 1041 t = t->_right; 1042 } 1043 else break; 1044 } 1045 1046 _leftTreeMax->_right = t->_left; 1047 _rightTreeMin->_left = t->_right; 1048 t->_left = header._right; 1049 t->_right = header._left; 1050 1051 } 1052 1053 //-------------------------------------------------------------------- 1054 //Rotate binary tree node with _left child. 1055 //-------------------------------------------------------------------- 1056 template <class T, class KeyType> 1057 void SplayTree<T, KeyType>::rotateWithLeftChild( BTreeNode<T, KeyType> * & k2 ) const 1058 { 1059 BTreeNode<T, KeyType> *k1 = k2->_left; 1060 k2->_left = k1->_right; 1061 k1->_right = k2; 1062 k2 = k1; 1063 } 1064 1065 //--------------------------------------------------------------------- 1066 //Rotate binary tree node with _right child. 1067 //--------------------------------------------------------------------- 1068 template <class T, class KeyType> 1069 void SplayTree<T, KeyType>::rotateWithRightChild( BTreeNode<T, KeyType> * & k1 ) const 1070 { 1071 BTreeNode<T, KeyType> *k2 = k1->_right; 1072 k1->_right = k2->_left; 1073 k2->_left = k1; 1074 k1 = k2; 1075 } 1076 1077 //---------------------------------------------------------------------- 1078 //Internal method to reclaim internal nodes in subtree t. 1079 //WARNING: This is prone to running out of stack space. 1080 //---------------------------------------------------------------------- 1081 template <class T, class KeyType> 1082 void SplayTree<T, KeyType>::reclaimMemory( BTreeNode<T, KeyType> * t ) const 1083 { 1084 if( t != t->_left ) 1085 { 1086 reclaimMemory( t->_left ); 1087 reclaimMemory( t->_right ); 1088 delete t; 1089 } 1090 } 1091 1092 //----------------------------------------------------------------------- 1093 //Internal method to clone subtree. 1094 //WARNING: This is prone to running out of stack space. 1095 //----------------------------------------------------------------------- 1096 template <class T, class KeyType> 1097 BTreeNode<T, KeyType> * SplayTree<T, KeyType>::clone( BTreeNode<T, KeyType> * t ) const 1098 { 1099 if( t == t->_left ) // Cannot test against NULLNode!!! 1100 return NULL; 1101 else 1102 return new BTreeNode<T, KeyType>( t->_data, clone( t->_left ), clone( t->_right ) ); 1103 } 1104 1105 //----------------------------------------------------------------------- 1106 //Tranverse the tree by pre-order method; 1107 //----------------------------------------------------------------------- 1108 template<class T, class KeyType> 1109 void SplayTree<T, KeyType>::PreOrder( void(*Visit)(BTreeNode<T, KeyType> *u), BTreeNode<T, KeyType> *t) 1110 { 1111 if(t!=NULL) 1112 { 1113 Visit(t); 1114 PreOrder(Visit,t->_left); 1115 PreOrder(Visit,t->_right); 1116 } 1117 1118 } 1119 1120 //----------------------------------------------------------------------- 1121 //Tranverse the tree by in-order method; 1122 //----------------------------------------------------------------------- 1123 template<class T, class KeyType> 1124 void SplayTree<T, KeyType>::InOrder( void(*Visit)(BTreeNode<T, KeyType> *u), BTreeNode<T, KeyType> *t) 1125 { 1126 if(t!=NULL) 1127 { 1128 InOrder(Visit,t->_left); 1129 Visit(t); 1130 InOrder(Visit,t->_right); 1131 } 1132 } 1133 1134 1135 //----------------------------------------------------------------------- 1136 //Tranverse the tree by in-order method; 1137 //----------------------------------------------------------------------- 1138 template<class T, class KeyType> 1139 void SplayTree<T, KeyType>::InOrder( void(*Visit)(BTreeNode<T, KeyType>*u, double y) 1140 , BTreeNode<T, KeyType> *t, double y) 1141 { 1142 if(t!=NULL) 1143 { 1144 InOrder(Visit,t->_left, y); 1145 Visit(t, y); 1146 InOrder(Visit,t->_right, y); 1147 } 1148 } 1149 1150 1151 1152 //----------------------------------------------------------------------- 1153 //Tranverse the tree by post-order method; 1154 //----------------------------------------------------------------------- 1155 template<class T, class KeyType> 1156 void SplayTree<T, KeyType>::PostOrder( void(*Visit)(BTreeNode<T, KeyType> *u), BTreeNode<T, KeyType> *t) 1157 { 1158 if(t!=NULL) 1159 { 1160 PostOrder(Visit,t->_left); 1161 PostOrder(Visit,t->_right); 1162 Visit(t); 1163 } 1164 } 1165 1166 //----------------------------------------------------------------------- 1167 //return the height of subtree 1168 //----------------------------------------------------------------------- 1169 template<class T, class KeyType> 1170 int SplayTree<T, KeyType>::Height(BTreeNode<T, KeyType> *subtree) const 1171 { 1172 if(subtree==NULL) return 0; 1173 int lh=Height(subtree->_left); 1174 int rh=Height(subtree->_right); 1175 1176 return (lh>rh)?(++lh):(++rh); 1177 } 1178 1179 1180 //////////////////////////////////////////////////////////////////////////////////////////// 1181 //////////////////////////////////////////////////////////////////////////////////////////// 1182 ////////////////////////////////// geometry.h ////////////////////////////////////////////// 1183 //////////////////////////////////////////////////////////////////////////////////////////// 1184 //////////////////////////////////////////////////////////////////////////////////////////// 1185 1186 //-------------------------------------------------------------------------/ 1187 //Copyright (C) 2003, 2004, 2005, ALL RIGHTS RESERVED. 1188 //Centre for Sys. Eng. & App. Mech. FEMAGSoft S.A. 1189 //Universite Cathalique de Louvain 4, Avenue Albert Einstein 1190 //Batiment Euler, Avenue Georges Lemaitre, 4 B-1348 Louvain-la-Neuve 1191 //B-1348, Louvain-la-Neuve Belgium 1192 //Belgium 1193 //-------------------------------------------------------------------------/ 1194 // 1195 //Name: geometry.h (all geometry premitives related polygon triang- 1196 // ulation by sweep line algorithm) 1197 //Author: Liang, Wu (wu@mema.ucl.ac.be, wuliang@femagsoft.com) 1198 //Created: 03/2001 1199 //Modified: 10/2005. Modified and simplified only for polygon triangul- 1200 // ation purpose. 1201 //-------------------------------------------------------------------------/ 1202 1203 1204 //base class for points; 1205 class Pointbase 1206 { 1207 public: 1208 //constructors and destructor 1209 Pointbase() {} 1210 Pointbase(const Pointbase& pb); 1211 1212 Pointbase(const double& xx, const double& yy) 1213 :id(0), x(xx), y(yy), type(UNKNOWN) { } 1214 1215 Pointbase(const int& idd, const double& xx, const double& yy) 1216 :id(idd), x(xx), y(yy), type(UNKNOWN) { } 1217 1218 Pointbase(const double& xx, const double& yy, const Type& ttype) 1219 :id(0), x(xx), y(yy), type(ttype) { } 1220 1221 Pointbase(const int& idd, const double& xx, const double& yy, const Type& ttype) 1222 :id(idd),x(xx), y(yy), type(ttype) { } 1223 1224 //operator overloading 1225 friend bool operator==(const Pointbase&, const Pointbase&); 1226 friend bool operator>(const Pointbase&, const Pointbase&); 1227 friend bool operator<(const Pointbase&, const Pointbase&); 1228 friend bool operator!=(const Pointbase&, const Pointbase&); 1229 1230 //public data 1231 unsigned int id; //id of point; 1232 double x, y; //coordinates; 1233 Type type; //type of points; 1234 bool left; //left chain or not; 1235 1236 }; 1237 1238 1239 //base class for polygon boundary 1240 //Linebase class is a directed line segment with start/end point 1241 class Linebase 1242 { 1243 public: 1244 //constructors and destructor 1245 Linebase(); 1246 Linebase(Pointbase* ep1, Pointbase* ep2, const Type& type,long int & l_id); 1247 Linebase(const Linebase& line); 1248 ~Linebase() {}; 1249 1250 unsigned int id() const { return _id; } 1251 1252 //two end points 1253 Pointbase* endPoint(const int& i) const { return _endp[i]; } 1254 Type type() const { return _type; } 1255 double keyValue() const { return _key; } 1256 void setKeyValue(const double& y); 1257 //slightly increased the key to avoid duplicated key for searching tree. 1258 void increaseKeyValue(const double& diff) { _key+=diff; } 1259 //reverse a directed line segment; reversable only for inserted diagonals 1260 void reverse(); 1261 1262 //set and return helper of a directed line segment; 1263 void setHelper(const unsigned int& i) { _helper=i; } 1264 unsigned int helper() { return _helper; } 1265 1266 protected: 1267 unsigned int _id; //id of a line segment; 1268 Pointbase* _endp[2]; //two end points; 1269 1270 Type _type; //type of a line segement, input/insert 1271 double _key; //key of a line segment for splay tree searching 1272 unsigned int _helper; //helper of a line segemnt 1273 }; 1274 1275 //////////////////////////////////////////////////////////////////////////////////////////// 1276 //////////////////////////////////////////////////////////////////////////////////////////// 1277 ////////////////////////////////// geometry.cc ////////////////////////////////////////////// 1278 //////////////////////////////////////////////////////////////////////////////////////////// 1279 //////////////////////////////////////////////////////////////////////////////////////////// 1280 1281 //-------------------------------------------------------------------------/ 1282 //Copyright (C) 2003, 2004, 2005, ALL RIGHTS RESERVED. 1283 //Centre for Sys. Eng. & App. Mech. FEMAGSoft S.A. 1284 //Universite Cathalique de Louvain 4, Avenue Albert Einstein 1285 //Batiment Euler, Avenue Georges Lemaitre, 4 B-1348 Louvain-la-Neuve 1286 //B-1348, Louvain-la-Neuve Belgium 1287 //Belgium 1288 //-------------------------------------------------------------------------/ 1289 // 1290 //Name: geometry.cc (all geometry premitive implementations related 1291 // to polygon triangulation by sweep line algorithm) 1292 //Author: Liang, Wu (wu@mema.ucl.ac.be, wuliang@femagsoft.com) 1293 //Created: 03/2001 1294 //Modified: 10/2005. Modified and simplified only for polygon triangul- 1295 // ation purpose. 1296 //-------------------------------------------------------------------------/ 1297 1298 1299 //Jonathan schewchuk's exact arithmetic code, see predicates.cc for detais; 1300 extern double orient2d(double* pa, double* pb, double* pc); 1301 1302 //---------------------------------------------------------------------------- 1303 //square of the distance of two points; 1304 //---------------------------------------------------------------------------- 1305 double dist_sqr(const Pointbase& sp, const Pointbase& ep) 1306 { 1307 return sqr(sp.x-ep.x)+sqr(sp.y-ep.y); 1308 } 1309 1310 //---------------------------------------------------------------------------- 1311 //square of the distance of two points; 1312 //---------------------------------------------------------------------------- 1313 double dist_sqr(double *pa, double *pb) 1314 { 1315 return sqr(pa[0]-pb[0])+sqr(pa[1]-pb[1]); 1316 } 1317 1318 void UpdateKey(BTreeNode<Linebase*,double>* node, double y) 1319 { 1320 node->data()->setKeyValue(y); 1321 } 1322 1323 //---------------------------------------------------------------------------- 1324 //copy constructor 1325 //---------------------------------------------------------------------------- 1326 Pointbase::Pointbase(const Pointbase& pb) 1327 { 1328 this->id=pb.id; 1329 this->x=pb.x; 1330 this->y=pb.y; 1331 this->type=pb.type; 1332 this->left=pb.left; 1333 } 1334 1335 //---------------------------------------------------------------------------- 1336 //operator ( ==, >, < and != ) overloading for pointbase class 1337 //---------------------------------------------------------------------------- 1338 bool operator==(const Pointbase& pa, const Pointbase& pb) 1339 { 1340 return (pa.x==pb.x && pa.y==pb.y); 1341 } 1342 1343 //---------------------------------------------------------------------------- 1344 bool operator>(const Pointbase& pa, const Pointbase& pb) 1345 { 1346 return( (pa.y > pb.y) || ( (pa.y==pb.y) && (pa.x < pb.x)) ); 1347 } 1348 1349 //---------------------------------------------------------------------------- 1350 bool operator<(const Pointbase& pa, const Pointbase& pb) 1351 { 1352 return( (pa.y < pb.y) || ( (pa.y==pb.y) && (pa.x > pb.x)) ); 1353 } 1354 1355 //---------------------------------------------------------------------------- 1356 bool operator!=(const Pointbase& pa, const Pointbase& pb) 1357 { 1358 return !(pa.x==pb.x && pa.y==pb.y); 1359 } 1360 1361 //---------------------------------------------------------------------------- 1362 //Linebase construct 1363 //---------------------------------------------------------------------------- 1364 Linebase::Linebase():_type(UNKNOWN) 1365 { 1366 for(int i=0; i<2; ++i) _endp[i]=0; 1367 _id=0; 1368 } 1369 1370 //----------------------------------------------------------------------------- 1371 //Linebase construct 1372 //----------------------------------------------------------------------------- 1373 Linebase::Linebase(Pointbase* sp, Pointbase* ep, const Type& type, long int & l_id):_type(type) 1374 { 1375 _endp[0]=sp; 1376 _endp[1]=ep; 1377 _id=++l_id; 1378 } 1379 1380 //---------------------------------------------------------------------------- 1381 //copy constructor 1382 //---------------------------------------------------------------------------- 1383 Linebase::Linebase(const Linebase& line) 1384 { 1385 this->_id=line._id; 1386 this->_endp[0]=line._endp[0]; 1387 this->_endp[1]=line._endp[1]; 1388 this->_key=line._key; 1389 this->_helper=line._helper; 1390 } 1391 1392 1393 //---------------------------------------------------------------------------- 1394 //reverse a directed line segment, reverseable only for insert diagonals 1395 //---------------------------------------------------------------------------- 1396 void Linebase::reverse() 1397 { 1398 assert(_type==INSERT); 1399 Pointbase* tmp=_endp[0]; 1400 _endp[0]=_endp[1]; 1401 _endp[1]=tmp; 1402 } 1403 1404 void Linebase::setKeyValue(const double& y) 1405 { 1406 if( _endp[1]->y==_endp[0]->y ) 1407 _key=_endp[0]->x < _endp[1]->x ? _endp[0]->x:_endp[1]->x; 1408 else _key=( y - _endp[0]->y ) * ( _endp[1]->x - _endp[0]->x ) / (_endp[1]->y - _endp[0]->y ) + _endp[0]->x; 1409 } 1410 1411 }//end namespace internal_poltrig 1412 1413 class PolygonTriangulator::Polygon 1414 { 1415 public: 1416 //constructor and destructor 1417 Polygon(const std::vector<double>& x,const std::vector<double>& y); 1418 ~Polygon(); 1419 1420 // main member function for polygon triangulation; 1421 void partition2Monotone(); 1422 void searchMonotones(); 1423 void triangulation(); 1424 1425 //return all triangles 1426 const Triangles* triangles() { return &_triangles; } 1427 1428 internal_poltrig::PointbaseMap& points() { return _points; } 1429 internal_poltrig::LineMap& edges() { return _edges; } 1430 1431 //private member functions. 1432 private: 1433 long int l_id;//previous a global, but that gives mem. crashes. Must be reinitted for every polygon. 1434 1435 void set_contour (const std::vector<double>& x,const std::vector<double>& y); 1436 void initializate(); 1437 1438 //prev or next point/edge id for a given ith point/edge; 1439 unsigned int prev(const unsigned int& i); 1440 unsigned int next(const unsigned int& i); 1441 1442 //handle event vertext according to vertex type; 1443 void handleStartVertex(const unsigned int& ); 1444 void handleEndVertex(const unsigned int& ); 1445 void handleSplitVertex(const unsigned int& ); 1446 void handleMergeVertex(const unsigned int& ); 1447 void handleRegularVertexUp(const unsigned int& ); 1448 void handleRegularVertexDown(const unsigned int& ); 1449 1450 //add diagonal between two vertices; 1451 void addDiagonal(const unsigned int& i, const unsigned int& j); 1452 1453 1454 //angle ABC for three given points, for monotone polygon searching purpose; 1455 double angleCosb(double *A, double *B, double *C); 1456 //find the next edge, for monotone polygon searching purpose; 1457 unsigned int selectNextEdge(internal_poltrig::Linebase* edge); 1458 1459 //triangulate a monotone polygon piece; 1460 void triangulateMonotone(internal_poltrig::Monopoly& mpoly); 1461 1462 //private data memebers 1463 internal_poltrig::PQueue _qpoints; //priority queue for event points 1464 internal_poltrig::EdgeBST _edgebst; //edge binary searching tree (splaytree) 1465 internal_poltrig::Monopolys _mpolys; //all monotone polygon piece list; 1466 Triangles _triangles; //all triangle list; 1467 1468 //data for monotone piece searching purpose; 1469 internal_poltrig::AdjEdgeMap _startAdjEdgeMap; //all edges starting from given points (map) 1470 internal_poltrig::LineMap _diagonals; //added diagonals to partition polygon to 1471 //monotont pieces, not all diagonals of 1472 1473 void init_vertices_and_lines(); 1474 unsigned int _ncontours; //number of contours 1475 std::vector<unsigned int> _nVertices; // 1476 internal_poltrig::PointbaseMap _points; //all vertices 1477 internal_poltrig::LineMap _edges; //all edges 1478 double _xmin,_xmax, _ymin,_ymax; //boundary box for polygon 1479 1480 }; 1481 1482 1483 void PolygonTriangulator::Polygon::init_vertices_and_lines() 1484 { 1485 int sid,eid; 1486 int num=0; 1487 1488 internal_poltrig::Type type; 1489 1490 for(unsigned j=0; j<_ncontours; ++j) 1491 { 1492 for(unsigned i=1; i<=_nVertices[j]; ++i)//fixme: 0-based? 1493 { 1494 sid=num+i; 1495 eid=(i==_nVertices[j])?num+1:num+i+1; 1496 type = internal_poltrig::INPUT; 1497 internal_poltrig::Linebase* line=new internal_poltrig::Linebase(_points[sid], _points[eid], type,l_id); 1498 _edges[l_id]=line; 1499 } 1500 num+=_nVertices[j]; 1501 } 1502 1503 int sum=0; 1504 for(unsigned int i=0; i<_ncontours; ++i) 1505 { 1506 sum+= _nVertices[i]; 1507 _nVertices[i]=sum; 1508 } 1509 1510 } 1511 1512 1513 //---------------------------------------------------------------------------- 1514 //get input 1515 //---------------------------------------------------------------------------- 1516 void PolygonTriangulator::Polygon::set_contour(const std::vector<double>& x,const std::vector<double>& y) 1517 { 1518 assert(x.size()==y.size()); 1519 1520 _nVertices.push_back( x.size() ); 1521 unsigned int i = _points.size()+1/*1*/;//fixme: get rid of the +1 ? 1522 double xx,yy; 1523 internal_poltrig::Type type; 1524 for (unsigned int j = 0; j < _nVertices[_ncontours]; ++j, ++i) 1525 { 1526 xx=x.at(j); 1527 yy=y.at(j); 1528 type=internal_poltrig::INPUT; 1529 internal_poltrig::Pointbase* point=new internal_poltrig::Pointbase(i,xx,yy,type); 1530 if(xx > _xmax ) _xmax=xx; 1531 if(xx < _xmin ) _xmin=xx; 1532 if(yy > _ymax ) _ymax=yy; 1533 if(yy < _ymin ) _ymin=yy; 1534 _points[i]=point; 1535 } 1536 1537 ++_ncontours; 1538 1539 } 1540 1541 1542 1543 1544 //---------------------------------------------------------------------------- 1545 //polygon class constructor 1546 //---------------------------------------------------------------------------- 1547 PolygonTriangulator::Polygon::Polygon(const std::vector<double>& x,const std::vector<double>& y) 1548 { 1549 l_id = 0; 1550 _ncontours=0; 1551 1552 _xmin = _ymin = std::numeric_limits<double>::infinity(); 1553 _xmax = _ymax = - std::numeric_limits<double>::infinity(); 1554 1555 1556 set_contour(x,y); 1557 init_vertices_and_lines(); 1558 initializate(); 1559 } 1560 1561 //---------------------------------------------------------------------------- 1562 //polygon destructor 1563 //---------------------------------------------------------------------------- 1564 PolygonTriangulator::Polygon::~Polygon() 1565 { 1566 //clear all dynamic allocated memory 1567 internal_poltrig::PointbaseMap::iterator itp=_points.begin(); 1568 for(; itp!=_points.end(); ++itp) 1569 { 1570 delete itp->second; 1571 } 1572 1573 internal_poltrig::LineMap::iterator itl=_edges.begin(); 1574 for(; itl!=_edges.end(); ++itl) 1575 { 1576 delete itl->second; 1577 } 1578 1579 } 1580 1581 //---------------------------------------------------------------------------- 1582 //return the previous point (or edge) id for a given ith point (or edge); 1583 //---------------------------------------------------------------------------- 1584 unsigned int PolygonTriangulator::Polygon::prev(const unsigned int& i) 1585 { 1586 unsigned int j(0),prevLoop(0),currentLoop(0); 1587 1588 while ( i > _nVertices[currentLoop] ) 1589 { 1590 prevLoop=currentLoop; 1591 ++currentLoop; 1592 } 1593 1594 if( i==1 || (i==_nVertices[prevLoop]+1) ) j=_nVertices[currentLoop]; 1595 else if( i <= _nVertices[currentLoop] ) j=i-1; 1596 1597 return j; 1598 } 1599 1600 //---------------------------------------------------------------------------- 1601 //return the next point (or edge) id for a given ith point (or edge); 1602 //---------------------------------------------------------------------------- 1603 unsigned int PolygonTriangulator::Polygon::next(const unsigned int& i) 1604 { 1605 unsigned int j(0),prevLoop(0),currentLoop(0); 1606 1607 while ( i > _nVertices[currentLoop] ) 1608 { 1609 prevLoop=currentLoop; 1610 ++currentLoop; 1611 } 1612 1613 if( i < _nVertices[currentLoop] ) j=i+1; 1614 else if ( i==_nVertices[currentLoop] ) 1615 { 1616 if( currentLoop==0) j=1; 1617 else j=_nVertices[prevLoop]+1; 1618 } 1619 1620 return j; 1621 } 1622 1623 //---------------------------------------------------------------------------- 1624 //polygon initialization; 1625 //to find types of all polygon vertices; 1626 //create a priority queue for all vertices; 1627 //construct an edge set for each vertex (the set holds all edges starting from 1628 //the vertex, only for loop searching purpose). 1629 //---------------------------------------------------------------------------- 1630 void PolygonTriangulator::Polygon::initializate() 1631 { 1632 internal_poltrig::PointbaseMap::iterator it=_points.begin(); 1633 for(; it!=_points.end(); ++it) 1634 { 1635 int id=it->first; 1636 int idp=prev(id); 1637 int idn=next(id); 1638 internal_poltrig::Pointbase p=*_points[id], pnext=*_points[idn], pprev=*_points[idp]; 1639 1640 if( p > pnext && pprev > p ) 1641 _points[id]->type=internal_poltrig::REGULAR_DOWN; 1642 else if (p > pprev && pnext > p) 1643 _points[id]->type=internal_poltrig::REGULAR_UP; 1644 else 1645 { 1646 double pa[2], pb[2], pc[2]; 1647 1648 pa[0]=_points[idp]->x; 1649 pa[1]=_points[idp]->y; 1650 1651 pb[0]=_points[id]->x; 1652 pb[1]=_points[id]->y; 1653 1654 pc[0]=_points[idn]->x; 1655 pc[1]=_points[idn]->y; 1656 1657 double area=internal_poltrig::orient2d(pa,pb,pc); 1658 1659 if( pprev > p && pnext > p ) _points[id]->type=(area >0) ? internal_poltrig::END: internal_poltrig::MERGE ; 1660 if( pprev < p && pnext < p ) _points[id]->type=(area >0) ? internal_poltrig::START : internal_poltrig::SPLIT; 1661 1662 } 1663 1664 _qpoints.push(*(it->second)); 1665 1666 _startAdjEdgeMap[id].insert(id); 1667 1668 } 1669 } 1670 1671 //---------------------------------------------------------------------------- 1672 //Add a diagonal from point id i to j 1673 //---------------------------------------------------------------------------- 1674 void PolygonTriangulator::Polygon::addDiagonal(const unsigned int& i, const unsigned int& j) 1675 { 1676 internal_poltrig::Type type=internal_poltrig::INSERT; 1677 internal_poltrig::Linebase* diag=new internal_poltrig::Linebase(_points[i], _points[j], type,l_id); 1678 _edges[diag->id()]=diag; 1679 1680 _startAdjEdgeMap[i].insert(diag->id()); 1681 _startAdjEdgeMap[j].insert(diag->id()); 1682 1683 _diagonals[diag->id()]=diag; 1684 1685 } 1686 1687 //---------------------------------------------------------------------------- 1688 //Handle start vertex 1689 //---------------------------------------------------------------------------- 1690 void PolygonTriangulator::Polygon::handleStartVertex(const unsigned int& i) 1691 { 1692 double y=_points[i]->y; 1693 _edgebst.InOrder(internal_poltrig::UpdateKey, y); 1694 1695 _edges[i]->setHelper(i); 1696 _edges[i]->setKeyValue(y); 1697 _edgebst.Insert(_edges[i]); 1698 1699 } 1700 1701 //---------------------------------------------------------------------------- 1702 //Handle end vertex 1703 //---------------------------------------------------------------------------- 1704 void PolygonTriangulator::Polygon::handleEndVertex(const unsigned int& i) 1705 { 1706 double y=_points[i]->y; 1707 _edgebst.InOrder(internal_poltrig::UpdateKey, y); 1708 1709 unsigned int previ=prev(i); 1710 internal_poltrig::Linebase* edge=_edges[previ]; 1711 unsigned int helper=_edges[previ]->helper(); 1712 1713 1714 if(_points[helper]->type==internal_poltrig::MERGE) addDiagonal(i, helper); 1715 _edgebst.Delete(edge->keyValue()); 1716 1717 } 1718 1719 //---------------------------------------------------------------------------- 1720 //Handle split vertex 1721 //---------------------------------------------------------------------------- 1722 void PolygonTriangulator::Polygon::handleSplitVertex(const unsigned int& i) 1723 { 1724 double x=_points[i]->x, y=_points[i]->y; 1725 _edgebst.InOrder(internal_poltrig::UpdateKey, y); 1726 1727 internal_poltrig::BTreeNode<internal_poltrig::Linebase*, double>* leftnode; 1728 _edgebst.FindMaxSmallerThan(x, leftnode); 1729 internal_poltrig::Linebase* leftedge=leftnode->data(); 1730 1731 unsigned int helper=leftedge->helper(); 1732 addDiagonal(i, helper); 1733 1734 leftedge->setHelper(i); 1735 _edges[i]->setHelper(i); 1736 _edges[i]->setKeyValue(y); 1737 _edgebst.Insert(_edges[i]); 1738 } 1739 1740 //---------------------------------------------------------------------------- 1741 //Handle merge vertex 1742 //---------------------------------------------------------------------------- 1743 void PolygonTriangulator::Polygon::handleMergeVertex(const unsigned int& i) 1744 { 1745 double x=_points[i]->x, y=_points[i]->y; 1746 _edgebst.InOrder(internal_poltrig::UpdateKey, y); 1747 1748 unsigned int previ=prev(i); 1749 unsigned int helper=_edges[previ]->helper(); 1750 if (_points[helper]->type==internal_poltrig::MERGE) addDiagonal(i, helper); 1751 _edgebst.Delete(_edges[previ]->keyValue()); 1752 1753 internal_poltrig::BTreeNode<internal_poltrig::Linebase*, double>* leftnode; 1754 _edgebst.FindMaxSmallerThan(x, leftnode); 1755 internal_poltrig::Linebase* leftedge=leftnode->data(); 1756 1757 helper=leftedge->helper(); 1758 if(_points[helper]->type==internal_poltrig::MERGE) addDiagonal(i, helper); 1759 1760 leftedge->setHelper(i); 1761 } 1762 1763 //---------------------------------------------------------------------------- 1764 //Handle regular down vertex 1765 //---------------------------------------------------------------------------- 1766 void PolygonTriangulator::Polygon::handleRegularVertexDown(const unsigned int& i) 1767 { 1768 double y=_points[i]->y; 1769 _edgebst.InOrder(internal_poltrig::UpdateKey, y); 1770 1771 unsigned int previ=prev(i); 1772 unsigned int helper=_edges[previ]->helper(); 1773 if(_points[helper]->type==internal_poltrig::MERGE) addDiagonal(i, helper); 1774 1775 _edgebst.Delete(_edges[previ]->keyValue()); 1776 _edges[i]->setHelper(i); 1777 _edges[i]->setKeyValue(y); 1778 _edgebst.Insert(_edges[i]); 1779 1780 } 1781 1782 //---------------------------------------------------------------------------- 1783 ////Handle regular up vertex 1784 //---------------------------------------------------------------------------- 1785 void PolygonTriangulator::Polygon::handleRegularVertexUp(const unsigned int& i) 1786 { 1787 double x=_points[i]->x, y=_points[i]->y; 1788 _edgebst.InOrder(internal_poltrig::UpdateKey, y); 1789 1790 internal_poltrig::BTreeNode<internal_poltrig::Linebase*, double>* leftnode; 1791 _edgebst.FindMaxSmallerThan(x, leftnode); 1792 1793 internal_poltrig::Linebase* leftedge=leftnode->data(); 1794 1795 unsigned int helper=leftedge->helper(); 1796 if(_points[helper]->type==internal_poltrig::MERGE) addDiagonal(i, helper); 1797 leftedge->setHelper(i); 1798 1799 } 1800 1801 //---------------------------------------------------------------------------- 1802 //partition polygon to monotone polygon pieces 1803 //---------------------------------------------------------------------------- 1804 void PolygonTriangulator::Polygon::partition2Monotone() 1805 { 1806 1807 if(_qpoints.top().type!=internal_poltrig::START) 1808 { 1809 std::cout<<"Please check your input polygon:\n1)orientations?\n2)duplicated points?\n"; 1810 exit(1); 1811 } 1812 1813 while(!_qpoints.empty()) 1814 { 1815 internal_poltrig::Pointbase vertex=_qpoints.top(); 1816 _qpoints.pop(); 1817 unsigned int id=vertex.id; 1818 1819 switch(vertex.type) 1820 { 1821 case internal_poltrig::START: handleStartVertex(id); break; 1822 case internal_poltrig::END: handleEndVertex(id); break; 1823 case internal_poltrig::MERGE: handleMergeVertex(id); break; 1824 case internal_poltrig::SPLIT: handleSplitVertex(id); break; 1825 case internal_poltrig::REGULAR_UP: handleRegularVertexUp(id); break; 1826 case internal_poltrig::REGULAR_DOWN: handleRegularVertexDown(id); break; 1827 default: 1828 std::cout<<"No duplicated points please!\n"; 1829 exit(1); break; 1830 } 1831 } 1832 } 1833 1834 1835 //---------------------------------------------------------------------------- 1836 //two Auxiliary functions to find monotone polygon pieces 1837 //---------------------------------------------------------------------------- 1838 1839 //---------------------------------------------------------------------------- 1840 //calculate angle B for A, B, C three given points 1841 //---------------------------------------------------------------------------- 1842 double PolygonTriangulator::Polygon::angleCosb(double *pa, double *pb, double *pc) 1843 { 1844 double dxab = pa[0] - pb[0]; 1845 double dyab = pa[1] - pb[1]; 1846 1847 double dxcb = pc[0] - pb[0]; 1848 double dycb = pc[1] - pb[1]; 1849 1850 double dxab2 = dxab * dxab; 1851 double dyab2 = dyab * dyab; 1852 double dxcb2 = dxcb * dxcb; 1853 double dycb2 = dycb * dycb; 1854 double ab = dxab2 + dyab2; 1855 double cb = dxcb2 + dycb2; 1856 1857 double cosb = dxab * dxcb + dyab * dycb; 1858 double denom = sqrt( ab * cb); 1859 1860 cosb/=denom; 1861 1862 return cosb; 1863 } 1864 1865 //---------------------------------------------------------------------------- 1866 //for any given edge, find the next edge we should choose when searching for 1867 //monotone polygon pieces; 1868 //---------------------------------------------------------------------------- 1869 unsigned int PolygonTriangulator::Polygon::selectNextEdge(internal_poltrig::Linebase* edge) 1870 { 1871 1872 unsigned int eid= edge->endPoint(1)->id; 1873 std::set<unsigned int> edges=_startAdjEdgeMap[eid]; 1874 assert(!edges.empty()); 1875 1876 unsigned int nexte=0; 1877 if( edges.size() == 1 ) nexte=*(edges.begin()); 1878 else if( edges.size() > 1 ) 1879 { 1880 unsigned int nexte_ccw(0), nexte_cw(0); 1881 double max=-2.0,min=2.0; 1882 1883 1884 std::set<unsigned int>::iterator it=edges.begin(); 1885 for(; it!=edges.end(); ++it) 1886 { 1887 if(*it==edge->id()) continue; 1888 double A[2], B[2], C[2]; 1889 A[0]=edge->endPoint(0)->x; A[1]=edge->endPoint(0)->y; 1890 B[0]=edge->endPoint(1)->x; B[1]=edge->endPoint(1)->y; 1891 1892 if(edge->endPoint(1)!=_edges[*it]->endPoint(0)) _edges[*it]->reverse(); 1893 C[0]=_edges[*it]->endPoint(1)->x; C[1]=_edges[*it]->endPoint(1)->y; 1894 1895 double area=internal_poltrig::orient2d(A, B, C); 1896 double cosb=angleCosb(A, B, C); 1897 1898 if( area > 0 && max < cosb ) { nexte_ccw=*it; max=cosb; } 1899 else if( min > cosb ) { nexte_cw=*it; min=cosb; } 1900 } 1901 1902 nexte = (nexte_ccw!=0) ? nexte_ccw : nexte_cw; 1903 } 1904 1905 return nexte; 1906 } 1907 1908 //---------------------------------------------------------------------------- 1909 //searching all monotone pieces; 1910 //---------------------------------------------------------------------------- 1911 void PolygonTriangulator::Polygon::searchMonotones() 1912 { 1913 int loop=0; 1914 1915 internal_poltrig::LineMap edges=_edges; 1916 1917 while( edges.size() > _diagonals.size() ) 1918 { 1919 ++loop; 1920 internal_poltrig::Monopoly poly; 1921 internal_poltrig::LineMap::iterator it=edges.begin(); 1922 internal_poltrig::Pointbase* startp=it->second->endPoint(0); 1923 internal_poltrig::Pointbase* endp=0; 1924 internal_poltrig::Linebase* next=it->second; 1925 1926 poly.push_back(startp->id); 1927 1928 for(;;) 1929 { 1930 endp=next->endPoint(1); 1931 if(next->type()!=internal_poltrig::INSERT) 1932 { 1933 edges.erase(next->id()); 1934 _startAdjEdgeMap[next->