src/glu/sgi/libnurbs/nurbtess/partitionY.cc

   1 /*
   2 ** License Applicability. Except to the extent portions of this file are
   3 ** made subject to an alternative license as permitted in the SGI Free
   4 ** Software License B, Version 1.1 (the "License"), the contents of this
   5 ** file are subject only to the provisions of the License. You may not use
   6 ** this file except in compliance with the License. You may obtain a copy
   7 ** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
   8 ** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
   9 **
  10 ** http://oss.sgi.com/projects/FreeB
  11 **
  12 ** Note that, as provided in the License, the Software is distributed on an
  13 ** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
  14 ** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
  15 ** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
  16 ** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
  17 **
  18 ** Original Code. The Original Code is: OpenGL Sample Implementation,
  19 ** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
  20 ** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
  21 ** Copyright in any portions created by third parties is as indicated
  22 ** elsewhere herein. All Rights Reserved.
  23 **
  24 ** Additional Notice Provisions: The application programming interfaces
  25 ** established by SGI in conjunction with the Original Code are The
  26 ** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
  27 ** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
  28 ** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
  29 ** Window System(R) (Version 1.3), released October 19, 1998. This software
  30 ** was created using the OpenGL(R) version 1.2.1 Sample Implementation
  31 ** published by SGI, but has not been independently verified as being
  32 ** compliant with the OpenGL(R) version 1.2.1 Specification.
  33 **
  34 */
  35 /*
  36 */
  37
  38 #include <stdlib.h>
  39 #include <stdio.h>
  40 #include <time.h>
  41
  42 #include "zlassert.h"
  43 #include "partitionY.h"
  44 #include "searchTree.h"
  45 #include "quicksort.h"
  46 #include "polyUtil.h"
  47
  48
  49 #define max(a,b) ((a>b)? a:b)
  50 #define min(a,b) ((a>b)? b:a)
  51
  52
  53 /*retrurn
  54  *-1: if A < B (Ya<Yb) || (Ya==Yb)
  55  * 0: if A == B
  56  * 1: if A>B
  57  */
  58 static Int compVertInY(Real A[2], Real B[2])
  59 {
  60   if( (A[1] < B[1]) || (A[1]==B[1] && A[0]<B[0]))
  61     return -1;
  62   else if
  63     ( A[1] == B[1] && A[0] == B[0]) return 0;
  64   else
  65     return 1;
  66 }
  67
  68 /*v is a vertex: the head of en edge,
  69  *e is an edge,
  70  *return 1 if e is below v: assume v1 and v2 are the two endpoints of e:
  71  * v1<= v, v2<=v.
  72  */
  73 Int isBelow(directedLine *v, directedLine *e)
  74 {
  75   Real* vert = v->head();
  76   if(   compVertInY(e->head(), vert) != 1
  77      && compVertInY(e->tail(), vert) != 1
  78      )
  79     return 1;
  80   else
  81     return 0;
  82 }
  83
  84 /*v is a vertex: the head of en edge,
  85  *e is an edge,
  86  *return 1 if e is below v: assume v1 and v2 are the two endpoints of e:
  87  * v1>= v, v2>=v.
  88  */
  89 Int isAbove(directedLine *v, directedLine *e)
  90 {
  91   Real* vert = v->head();
  92   if(   compVertInY(e->head(), vert) != -1
  93      && compVertInY(e->tail(), vert) != -1
  94      )
  95     return 1;
  96   else
  97     return 0;
  98 }
  99
 100 Int isCusp(directedLine *v)
 101 {
 102   Real *A=v->getPrev()->head();
 103   Real *B=v->head();
 104   Real *C=v->tail();
 105   if(A[1] < B[1] && B[1] < C[1])
 106     return 0;
 107   else if(A[1] > B[1] && B[1] > C[1])
 108     return 0;
 109   else if(A[1] < B[1] && C[1] < B[1])
 110     return 1;
 111   else if(A[1] > B[1] && C[1] > B[1])
 112     return 1;
 113
 114   if((isAbove(v, v) && isAbove(v, v->getPrev())) ||
 115      (isBelow(v, v) && isBelow(v, v->getPrev())))
 116     return 1;
 117   else
 118     return 0;
 119 }
 120
 121 /*crossproduct is strictly less than 0*/
 122 Int isReflex(directedLine *v)
 123 {
 124   Real* A = v->getPrev()->head();
 125   Real* B = v->head();
 126   Real* C = v->tail();
 127   Real Bx,By, Cx, Cy;
 128   Bx = B[0] - A[0];
 129   By = B[1] - A[1];
 130   Cx = C[0] - A[0];
 131   Cy = C[1] - A[1];
 132
 133   if(Bx*Cy - Cx*By < 0) return 1;
 134   else return 0;
 135 }
 136
 137  /*return
 138  *0: not-cusp
 139  *1: interior cusp
 140  *2: exterior cusp
 141  */
 142 Int cuspType(directedLine *v)
 143 {
 144   if(! isCusp(v)) return 0;
 145   else if(isReflex(v)) return 1;
 146   else
 147     return 2;
 148 }
 149
 150 sweepRange* sweepRangeMake(directedLine* left, Int leftType,
 151                            directedLine* right, Int rightType)
 152 {
 153   sweepRange* ret = (sweepRange*)malloc(sizeof(sweepRange));
 154   assert(ret);
 155   ret->left = left;
 156   ret->leftType = leftType;
 157   ret->right = right;
 158   ret->rightType = rightType;
 159   return ret;
 160 }
 161
 162 void sweepRangeDelete(sweepRange* range)
 163 {
 164   free(range);
 165 }
 166
 167 Int sweepRangeEqual(sweepRange* src1, sweepRange* src2)
 168 {
 169   Int leftEqual;
 170   Int rightEqual;
 171
 172
 173   /*The case when both are vertices should not happen*/
 174   assert(! (src1->leftType == 0 && src2->leftType == 0));
 175   if(src1->leftType == 0 && src2->leftType == 1){
 176     if(src1->left == src2->left ||
 177        src1->left->getPrev() == src2->left
 178        )
 179       leftEqual = 1;
 180     else
 181       leftEqual = 0;
 182   }
 183   else if(src1->leftType == 1 && src2->leftType == 1){
 184     if(src1->left == src2->left)
 185       leftEqual = 1;
 186     else
 187       leftEqual = 0;
 188   }
 189   else /*src1->leftType == 1 && src2->leftType == 0*/{
 190     if(src1->left == src2->left ||
 191        src1->left == src2->left->getPrev()
 192        )
 193       leftEqual = 1;
 194     else
 195       leftEqual = 0;
 196   }
 197
 198   /*the same thing for right*/
 199   /*The case when both are vertices should not happen*/
 200   assert(! (src1->rightType == 0 && src2->rightType == 0));
 201   if(src1->rightType == 0 && src2->rightType == 1){
 202     if(src1->right == src2->right ||
 203        src1->right->getPrev() == src2->right
 204        )
 205       rightEqual = 1;
 206     else
 207       rightEqual = 0;
 208   }
 209   else if(src1->rightType == 1 && src2->rightType == 1){
 210     if(src1->right == src2->right)
 211       rightEqual = 1;
 212     else
 213       rightEqual = 0;
 214   }
 215   else /*src1->rightType == 1 && src2->rightType == 0*/{
 216     if(src1->right == src2->right ||
 217        src1->right == src2->right->getPrev()
 218        )
 219       rightEqual = 1;
 220     else
 221       rightEqual = 0;
 222   }
 223
 224   return (leftEqual == 1 || rightEqual == 1);
 225 }
 226
 227 /*given (x_1, y_1) and (x_2, y_2), and y
 228  *return x such that (x,y) is on the line
 229  */
 230 inline/*static*/ Real intersectHoriz(Real x1, Real y1, Real x2, Real y2, Real y)
 231 {
 232   return ((y2==y1)? (x1+x2)*Real(0.5) : x1 + ((y-y1)/(y2-y1)) * (x2-x1));
 233 /*
 234   if(y2 == y1) return (x1+x2)*0.5;
 235   else return x1 + ((y-y1)/(y2-y1)) * (x2-x1);
 236 */
 237 }
 238
 239 /*compare two edges of a polygon.
 240  *edge A < edge B if there is a horizontal line so that the intersection
 241  *with A is to the left of the intersection with B.
 242  *This function is used in sweepY for the dynamic search tree insertion to
 243  *order the edges.
 244  * Implementation: (x_1,y_1) and (x_2, y_2)
 245  */
 246 static Int compEdges(directedLine *e1, directedLine *e2)
 247 {
 248   Real* head1 = e1->head();
 249   Real* tail1 = e1->tail();
 250   Real* head2 = e2->head();
 251   Real* tail2 = e2->tail();
 252 /*
 253   Real h10 = head1[0];
 254   Real h11 = head1[1];
 255   Real t10 = tail1[0];
 256   Real t11 = tail1[1];
 257   Real h20 = head2[0];
 258   Real h21 = head2[1];
 259   Real t20 = tail2[0];
 260   Real t21 = tail2[1];
 261 */
 262   Real e1_Ymax, e1_Ymin, e2_Ymax, e2_Ymin;
 263 /*
 264   if(h11>t11) {
 265     e1_Ymax= h11;
 266     e1_Ymin= t11;
 267   }
 268   else{
 269     e1_Ymax = t11;
 270     e1_Ymin = h11;
 271   }
 272
 273   if(h21>t21) {
 274     e2_Ymax= h21;
 275     e2_Ymin= t21;
 276   }
 277   else{
 278     e2_Ymax = t21;
 279     e2_Ymin = h21;
 280   }
 281 */
 282
 283   if(head1[1]>tail1[1]) {
 284     e1_Ymax= head1[1];
 285     e1_Ymin= tail1[1];
 286   }
 287   else{
 288     e1_Ymax = tail1[1];
 289     e1_Ymin = head1[1];
 290   }
 291
 292   if(head2[1]>tail2[1]) {
 293     e2_Ymax= head2[1];
 294     e2_Ymin= tail2[1];
 295   }
 296   else{
 297     e2_Ymax = tail2[1];
 298     e2_Ymin = head2[1];
 299   }
 300
 301
 302   /*Real e1_Ymax = max(head1[1], tail1[1]);*/ /*max(e1->head()[1], e1->tail()[1]);*/
 303   /*Real e1_Ymin = min(head1[1], tail1[1]);*/ /*min(e1->head()[1], e1->tail()[1]);*/
 304   /*Real e2_Ymax = max(head2[1], tail2[1]);*/ /*max(e2->head()[1], e2->tail()[1]);*/
 305   /*Real e2_Ymin = min(head2[1], tail2[1]);*/ /*min(e2->head()[1], e2->tail()[1]);*/
 306
 307   Real Ymax = min(e1_Ymax, e2_Ymax);
 308   Real Ymin = max(e1_Ymin, e2_Ymin);
 309
 310   Real y = Real(0.5)*(Ymax + Ymin);
 311
 312 /*  Real x1 = intersectHoriz(e1->head()[0], e1->head()[1], e1->tail()[0], e1->tail()[1], y);
 313   Real x2 = intersectHoriz(e2->head()[0], e2->head()[1], e2->tail()[0], e2->tail()[1], y);
 314 */
 315 /*
 316   Real x1 = intersectHoriz(h10, h11, t10, t11, y);
 317   Real x2 = intersectHoriz(h20, h21, t20, t21, y);
 318 */
 319   Real x1 = intersectHoriz(head1[0], head1[1], tail1[0], tail1[1], y);
 320   Real x2 = intersectHoriz(head2[0], head2[1], tail2[0], tail2[1], y);
 321
 322   if(x1<= x2) return -1;
 323   else return 1;
 324 }
 325
 326 /*used by sort precedures
 327  */
 328 static Int compInY(directedLine* v1, directedLine* v2)
 329 {
 330   return v1->compInY(v2);
 331 }
 332
 333 void findDiagonals(Int total_num_edges, directedLine** sortedVertices, sweepRange** ranges, Int& num_diagonals, directedLine** diagonal_vertices)
 334 {
 335   Int i,j,k;
 336
 337   k=0;
 338
 339   for(i=0; i<total_num_edges; i++)
 340     {
 341       directedLine* vert =sortedVertices[i];
 342       directedLine* thisEdge = vert;
 343       directedLine* prevEdge = vert->getPrev();
 344 /*
 345 printf("find i=%i\n", i);
 346 printf("the vertex is\n");
 347 vert->printSingle();
 348 */
 349       if(isBelow(vert, thisEdge) && isBelow(vert, prevEdge) && compEdges(prevEdge, thisEdge)<0)
 350         {
 351           /*this is an upward interior cusp*/
 352           diagonal_vertices[k++] = vert;
 353
 354           for(j=i+1; j<total_num_edges; j++)
 355             if(sweepRangeEqual(ranges[i], ranges[j]))
 356               {
 357                 diagonal_vertices[k++] = sortedVertices[j];
 358                 break;
 359               }
 360           assert(j<total_num_edges);
 361
 362
 363         }
 364       else if(isAbove(vert, thisEdge) && isAbove(vert, prevEdge) && compEdges(prevEdge, thisEdge)>0)
 365         {
 366           /*this is an downward interior cusp*/
 367           diagonal_vertices[k++] = vert;
 368           for(j=i-1; j>=0; j--)
 369             if(sweepRangeEqual(ranges[i], ranges[j]))
 370               {
 371                 diagonal_vertices[k++] = sortedVertices[j];
 372                 break;
 373               }
 374 /*        printf("j=%i\n", j);*/
 375           assert(j>=0);
 376
 377
 378
 379         }
 380     }
 381   num_diagonals = k/2;
 382 }
 383
 384 /*get rid of repeated diagonlas so that each diagonal appears only once in the array
 385  */
 386 Int deleteRepeatDiagonals(Int num_diagonals, directedLine** diagonal_vertices, directedLine** new_vertices)
 387 {
 388   Int i,k;
 389   Int j,l;
 390   Int index;
 391   index=0;
 392   for(i=0,k=0; i<num_diagonals; i++, k+=2)
 393     {
 394       Int isRepeated=0;
 395       /*check the diagonla (diagonal_vertice[k], diagonal_vertices[k+1])
 396        *is repeated or not
 397        */
 398       for(j=0,l=0; j<index; j++, l+=2)
 399         {
 400           if(
 401              (diagonal_vertices[k] == new_vertices[l] &&
 402               diagonal_vertices[k+1] == new_vertices[l+1]
 403               )
 404              ||
 405              (
 406               diagonal_vertices[k] == new_vertices[l+1] &&
 407               diagonal_vertices[k+1] == new_vertices[l]
 408               )
 409              )
 410             {
 411               isRepeated=1;
 412               break;
 413             }
 414         }
 415       if(! isRepeated)
 416         {
 417           new_vertices[index+index] = diagonal_vertices[k];
 418           new_vertices[index+index+1] = diagonal_vertices[k+1];
 419           index++;
 420         }
 421     }
 422   return index;
 423 }
 424
 425 /*for debug only*/
 426 directedLine** DBGfindDiagonals(directedLine *polygons, Int& num_diagonals)
 427 {
 428   Int total_num_edges = 0;
 429   directedLine** array = polygons->toArrayAllPolygons(total_num_edges);
 430   quicksort( (void**)array, 0, total_num_edges-1, (Int (*)(void*, void*)) compInY);
 431   sweepRange** ranges = (sweepRange**) malloc(sizeof(sweepRange*) * total_num_edges);
 432   assert(ranges);
 433
 434   sweepY(total_num_edges, array, ranges);
 435
 436  directedLine** diagonal_vertices = (directedLine**) malloc(sizeof(directedLine*) * total_num_edges);
 437   assert(diagonal_vertices);
 438   findDiagonals(total_num_edges, array, ranges, num_diagonals, diagonal_vertices);
 439
 440   num_diagonals=deleteRepeatDiagonals(num_diagonals, diagonal_vertices, diagonal_vertices);
 441   return diagonal_vertices;
 442
 443 }
 444
 445
 446 /*partition into Y-monotone polygons*/
 447 directedLine* partitionY(directedLine *polygons, sampledLine **retSampledLines)
 448 {
 449   Int total_num_edges = 0;
 450   directedLine** array = polygons->toArrayAllPolygons(total_num_edges);
 451
 452   quicksort( (void**)array, 0, total_num_edges-1, (Int (*)(void*, void*)) compInY);
 453
 454   sweepRange** ranges = (sweepRange**) malloc(sizeof(sweepRange*) * (total_num_edges));
 455   assert(ranges);
 456
 457
 458
 459   sweepY(total_num_edges, array, ranges);
 460
 461
 462
 463   /*the diagonal vertices are stored as:
 464    *v0-v1: 1st diagonal
 465    *v2-v3: 2nd diagonal
 466    *v5-v5: 3rd diagonal
 467    *...
 468    */
 469
 470
 471   Int num_diagonals;
 472   /*number diagonals is < total_num_edges*total_num_edges*/
 473   directedLine** diagonal_vertices = (directedLine**) malloc(sizeof(directedLine*) * total_num_edges*2/*total_num_edges*/);
 474   assert(diagonal_vertices);
 475
 476
 477
 478   findDiagonals(total_num_edges, array, ranges, num_diagonals, diagonal_vertices);
 479
 480
 481
 482   directedLine* ret_polygons = polygons;
 483   sampledLine* newSampledLines = NULL;
 484   Int i,k;
 485
 486 num_diagonals=deleteRepeatDiagonals(num_diagonals, diagonal_vertices, diagonal_vertices);
 487
 488
 489
 490   Int *removedDiagonals=(Int*)malloc(sizeof(Int) * num_diagonals);
 491   for(i=0; i<num_diagonals; i++)
 492     removedDiagonals[i] = 0;
 493
 494
 495
 496
 497
 498   for(i=0,k=0; i<num_diagonals; i++,k+=2)
 499     {
 500
 501
 502       directedLine* v1=diagonal_vertices[k];
 503       directedLine* v2=diagonal_vertices[k+1];
 504       directedLine* ret_p1;
 505       directedLine* ret_p2;
 506
 507       /*we ahve to determine whether v1 and v2 belong to the same polygon before
 508        *their structure are modified by connectDiagonal().
 509        */
 510 /*
 511       directedLine *root1 = v1->findRoot();
 512       directedLine *root2 = v2->findRoot();
 513       assert(root1);
 514       assert(root2);
 515 */
 516
 517 directedLine*  root1 = v1->rootLinkFindRoot();
 518 directedLine*  root2 = v2->rootLinkFindRoot();
 519
 520       if(root1 != root2)
 521         {
 522
 523           removedDiagonals[i] = 1;
 524           sampledLine* generatedLine;
 525
 526
 527
 528           v1->connectDiagonal(v1,v2, &ret_p1, &ret_p2, &generatedLine, ret_polygons);
 529
 530
 531
 532           newSampledLines = generatedLine->insert(newSampledLines);
 533 /*
 534           ret_polygons = ret_polygons->cutoffPolygon(root1);
 535
 536           ret_polygons = ret_polygons->cutoffPolygon(root2);
 537           ret_polygons = ret_p1->insertPolygon(ret_polygons);
 538 root1->rootLinkSet(ret_p1);
 539 root2->rootLinkSet(ret_p1);
 540 ret_p1->rootLinkSet(NULL);
 541 ret_p2->rootLinkSet(ret_p1);
 542 */
 543           ret_polygons = ret_polygons->cutoffPolygon(root2);
 544
 545
 546
 547 root2->rootLinkSet(root1);
 548 ret_p1->rootLinkSet(root1);
 549 ret_p2->rootLinkSet(root1);
 550
 551        /*now that we have connected the diagonal v1 and v2,
 552         *we have to check those unprocessed diagonals which
 553         *have v1 or v2 as an end point. Notice that the head of v1
 554         *has the same coodinates as the head of v2->prev, and the head of
 555         *v2 has the same coordinate as the head of v1->prev.
 556         *Suppose these is a diagonal (v1, x). If (v1,x) is still a valid
 557         *diagonal, then x should be on the left hand side of the directed line:        *v1->prev->head -- v1->head -- v1->tail. Otherwise, (v1,x) should be
 558         *replaced by (v2->prev, x), that is, x is on the left of
 559         * v2->prev->prev->head, v2->prev->head, v2->prev->tail.
 560         */
 561         Int ii, kk;
 562         for(ii=0, kk=0; ii<num_diagonals; ii++, kk+=2)
 563           if( removedDiagonals[ii]==0)
 564             {
 565               directedLine* d1=diagonal_vertices[kk];
 566               directedLine* d2=diagonal_vertices[kk+1];
 567               /*check d1, and replace diagonal_vertices[kk] if necessary*/
 568               if(d1 == v1) {
 569                 /*check if d2 is to left of v1->prev->head:v1->head:v1->tail*/
 570                 if(! pointLeft2Lines(v1->getPrev()->head(),
 571                                      v1->head(), v1->tail(), d2->head()))
 572                   {
 573 /*
 574                     assert(pointLeft2Lines(v2->getPrev()->getPrev()->head(),
 575                                            v2->getPrev()->head(),
 576                                            v2->getPrev()->tail(), d2->head()));
 577 */
 578                     diagonal_vertices[kk] = v2->getPrev();
 579                   }
 580               }
 581               if(d1 == v2) {
 582                 /*check if d2 is to left of v2->prev->head:v2->head:v2->tail*/
 583                 if(! pointLeft2Lines(v2->getPrev()->head(),
 584                                      v2->head(), v2->tail(), d2->head()))
 585                   {
 586 /*
 587                     assert(pointLeft2Lines(v1->getPrev()->getPrev()->head(),
 588                                            v1->getPrev()->head(),
 589                                            v1->getPrev()->tail(), d2->head()));
 590 */
 591                     diagonal_vertices[kk] = v1->getPrev();
 592                   }
 593               }
 594               /*check d2 and replace diagonal_vertices[k+1] if necessary*/
 595               if(d2 == v1) {
 596                 /*check if d1 is to left of v1->prev->head:v1->head:v1->tail*/
 597                 if(! pointLeft2Lines(v1->getPrev()->head(),
 598                                      v1->head(), v1->tail(), d1->head()))
 599                   {
 600 /*                  assert(pointLeft2Lines(v2->getPrev()->getPrev()->head(),
 601                                            v2->getPrev()->head(),
 602                                            v2->getPrev()->tail(), d1->head()));
 603 */
 604                     diagonal_vertices[kk+1] = v2->getPrev();
 605                   }
 606               }
 607               if(d2 == v2) {
 608                 /*check if d1 is to left of v2->prev->head:v2->head:v2->tail*/
 609                 if(! pointLeft2Lines(v2->getPrev()->head(),
 610                                      v2->head(), v2->tail(), d1->head()))
 611                   {
 612 /*                  assert(pointLeft2Lines(v1->getPrev()->getPrev()->head(),
 613                                            v1->getPrev()->head(),
 614                                            v1->getPrev()->tail(), d1->head()));
 615 */
 616                     diagonal_vertices[kk+1] = v1->getPrev();
 617                   }
 618               }
 619             }
 620 }/*end if (root1 not equal to root 2)*/
 621 }
 622
 623   /*second pass,  now all diagoals should belong to the same polygon*/
 624
 625
 626
 627   for(i=0,k=0; i<num_diagonals; i++, k += 2)
 628     if(removedDiagonals[i] == 0)
 629       {
 630
 631
 632         directedLine* v1=diagonal_vertices[k];
 633         directedLine* v2=diagonal_vertices[k+1];
 634
 635
 636
 637         directedLine* ret_p1;
 638         directedLine* ret_p2;
 639
 640         /*we ahve to determine whether v1 and v2 belong to the same polygon before
 641          *their structure are modified by connectDiagonal().
 642          */
 643         directedLine *root1 = v1->findRoot();
 644 /*
 645         directedLine *root2 = v2->findRoot();
 646
 647
 648
 649         assert(root1);
 650         assert(root2);
 651         assert(root1 == root2);
 652   */
 653         sampledLine* generatedLine;
 654
 655
 656
 657         v1->connectDiagonal(v1,v2, &ret_p1, &ret_p2, &generatedLine, ret_polygons);
 658         newSampledLines = generatedLine->insert(newSampledLines);
 659
 660         ret_polygons = ret_polygons->cutoffPolygon(root1);
 661
 662         ret_polygons = ret_p1->insertPolygon(ret_polygons);
 663
 664         ret_polygons = ret_p2->insertPolygon(ret_polygons);
 665
 666
 667
 668         for(Int j=i+1; j<num_diagonals; j++)
 669           {
 670             if(removedDiagonals[j] ==0)
 671               {
 672
 673                 directedLine* temp1=diagonal_vertices[2*j];
 674                 directedLine* temp2=diagonal_vertices[2*j+1];
 675                if(temp1==v1 || temp1==v2 || temp2==v1 || temp2==v2)
 676                 if(! temp1->samePolygon(temp1, temp2))
 677                   {
 678                     /*if temp1 and temp2 are in different polygons,
 679                      *then one of them must be v1 or v2.
 680                      */
 681
 682
 683
 684                     assert(temp1==v1 || temp1 == v2 || temp2==v1 || temp2 ==v2);
 685                     if(temp1==v1)
 686                       {
 687                         diagonal_vertices[2*j] = v2->getPrev();
 688                       }
 689                     if(temp2==v1)
 690                       {
 691                         diagonal_vertices[2*j+1] = v2->getPrev();
 692                       }
 693                     if(temp1==v2)
 694                       {
 695                         diagonal_vertices[2*j] = v1->getPrev();
 696                       }
 697                     if(temp2==v2)
 698                       {
 699                         diagonal_vertices[2*j+1] = v1->getPrev();
 700                       }
 701                   }
 702               }
 703           }
 704
 705       }
 706
 707   /*clean up spaces*/
 708   free(array);
 709   free(ranges);
 710   free(diagonal_vertices);
 711   free(removedDiagonals);
 712
 713   *retSampledLines = newSampledLines;
 714   return ret_polygons;
 715 }
 716
 717 /*given a set of simple polygons where the interior
 718  *is decided by left-hand principle,
 719  *return a range (sight) for each vertex. This is called
 720  *Trapezoidalization.
 721  */
 722 void sweepY(Int nVertices, directedLine** sortedVertices, sweepRange** ret_ranges)
 723 {
 724   Int i;
 725   /*for each vertex in the sorted list, update the binary search tree.
 726    *and store the range information for each vertex.
 727    */
 728   treeNode* searchTree = NULL;
 729   for(i=0; i<nVertices;i++)
 730     {
 731
 732       directedLine* vert = sortedVertices[i];
 733
 734       directedLine* thisEdge = vert;
 735       directedLine* prevEdge = vert->getPrev();
 736
 737       if(isBelow(vert, thisEdge) && isAbove(vert, prevEdge))
 738         {
 739
 740           /*case 1: this < v < prev
 741            *the polygon is going down at v, the interior is to
 742            *the right hand side.
 743            * find the edge to the right of thisEdge for right range.
 744            * delete thisEdge
 745            * insert prevEdge
 746            */
 747           treeNode* thisNode = TreeNodeFind(searchTree, thisEdge, ( Int (*) (void *, void *))compEdges);
 748           assert(thisNode);
 749
 750           treeNode* succ = TreeNodeSuccessor(thisNode);
 751           assert(succ);
 752           searchTree = TreeNodeDeleteSingleNode(searchTree, thisNode);
 753           searchTree = TreeNodeInsert(searchTree, TreeNodeMake(prevEdge), ( Int (*) (void *, void *))compEdges);
 754
 755
 756           ret_ranges[i] = sweepRangeMake(vert, 0, (directedLine*) (succ->key), 1);
 757
 758         }
 759       else if(isAbove(vert, thisEdge) && isBelow(vert, prevEdge))
 760         {
 761
 762           /*case 2: this > v > prev
 763            *the polygon is going up at v, the interior is to
 764            *the left hand side.
 765            * find the edge to the left of thisEdge for left range.
 766            * delete prevEdge
 767            * insert thisEdge
 768            */
 769           treeNode* prevNode = TreeNodeFind(searchTree, prevEdge, ( Int (*) (void *, void *))compEdges);
 770           assert(prevNode);
 771           treeNode* pred = TreeNodePredecessor(prevNode);
 772           searchTree = TreeNodeDeleteSingleNode(searchTree, prevNode);
 773           searchTree = TreeNodeInsert(searchTree, TreeNodeMake(thisEdge), ( Int (*) (void *, void *))compEdges);
 774           ret_ranges[i] = sweepRangeMake((directedLine*)(pred->key), 1, vert, 0);
 775         }
 776       else if(isAbove(vert, thisEdge) && isAbove(vert, prevEdge))
 777         {
 778
 779           /*case 3: insert both edges*/
 780           treeNode* thisNode = TreeNodeMake(thisEdge);
 781           treeNode* prevNode = TreeNodeMake(prevEdge);
 782           searchTree = TreeNodeInsert(searchTree, thisNode, ( Int (*) (void *, void *))compEdges);
 783           searchTree = TreeNodeInsert(searchTree, prevNode, ( Int (*) (void *, void *))compEdges);
 784           if(compEdges(thisEdge, prevEdge)<0)  /*interior cusp*/
 785             {
 786
 787               treeNode* leftEdge = TreeNodePredecessor(thisNode);
 788               treeNode* rightEdge = TreeNodeSuccessor(prevNode);
 789               ret_ranges[i] = sweepRangeMake( (directedLine*) leftEdge->key, 1,
 790                                              (directedLine*) rightEdge->key, 1
 791                                              );
 792             }
 793           else /*exterior cusp*/
 794             {
 795
 796               ret_ranges[i] = sweepRangeMake( prevEdge, 1, thisEdge, 1);
 797             }
 798         }
 799       else if(isBelow(vert, thisEdge) && isBelow(vert, prevEdge))
 800         {
 801
 802           /*case 4: delete both edges*/
 803           treeNode* thisNode = TreeNodeFind(searchTree, thisEdge, ( Int (*) (void *, void *))compEdges);
 804           treeNode* prevNode = TreeNodeFind(searchTree, prevEdge, ( Int (*) (void *, void *))compEdges);
 805           if(compEdges(thisEdge, prevEdge)>0) /*interior cusp*/
 806             {
 807               treeNode* leftEdge = TreeNodePredecessor(prevNode);
 808               treeNode* rightEdge = TreeNodeSuccessor(thisNode);
 809               ret_ranges[i] = sweepRangeMake( (directedLine*) leftEdge->key, 1,
 810                                              (directedLine*) rightEdge->key, 1
 811                                              );
 812             }
 813           else /*exterior cusp*/
 814             {
 815               ret_ranges[i] = sweepRangeMake( thisEdge, 1, prevEdge, 1);
 816             }
 817           searchTree = TreeNodeDeleteSingleNode(searchTree, thisNode);
 818           searchTree = TreeNodeDeleteSingleNode(searchTree, prevNode);
 819         }
 820       else
 821         {
 822           fprintf(stderr,"error in partitionY.C, invalid case\n");
 823           printf("vert is\n");
 824           vert->printSingle();
 825           printf("thisEdge is\n");
 826           thisEdge->printSingle();
 827           printf("prevEdge is\n");
 828           prevEdge->printSingle();
 829
 830           exit(1);
 831         }
 832     }
 833
 834   /*finaly clean up space: delete the search tree*/
 835   TreeNodeDeleteWholeTree(searchTree);
 836 }