src/glu/sgi/libnurbs/nurbtess/monoChain.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 "gluos.h"
  39 #include <stdlib.h>
  40 #include <stdio.h>
  41 #include <GL/gl.h>
  42
  43 #include "glimports.h"
  44 #include "zlassert.h"
  45
  46 #include "monoChain.h"
  47 #include "quicksort.h"
  48 #include "searchTree.h"
  49 #include "polyUtil.h"
  50
  51 #ifndef max
  52 #define max(a,b) ((a>b)? a:b)
  53 #endif
  54 #ifndef min
  55 #define min(a,b) ((a>b)? b:a)
  56 #endif
  57
  58 extern Int isCusp(directedLine *v);
  59 extern Int deleteRepeatDiagonals(Int num_diagonals, directedLine** diagonal_vertices, directedLine** new_vertices);
  60
  61 //for debug purpose only
  62 #if 0 // UNUSED
  63 static void drawDiagonals(Int num_diagonals, directedLine** diagonal_vertices)
  64 {
  65   Int i;
  66   for(i=0; i<num_diagonals; i++)
  67     {
  68       glBegin(GL_LINE);
  69       glVertex2fv(diagonal_vertices[2*i]->head());
  70       glVertex2fv(diagonal_vertices[2*i+1]->head());
  71       glEnd();
  72     }
  73 }
  74 #endif
  75
  76 /*given (x_1, y_1) and (x_2, y_2), and y
  77  *return x such that (x,y) is on the line
  78  */
  79 inline Real intersectHoriz(Real x1, Real y1, Real x2, Real y2, Real y)
  80 {
  81   return ((y2==y1)? (x1+x2)*0.5 : x1 + ((y-y1)/(y2-y1)) * (x2-x1));
  82 }
  83
  84 //compare the heads of the two chains
  85 static int compChainHeadInY(monoChain* mc1, monoChain* mc2)
  86 {
  87   return compV2InY(mc1->getHead()->head(), mc2->getHead()->head());
  88 }
  89
  90 monoChain::monoChain(directedLine* cHead, directedLine* cTail)
  91 {
  92   chainHead = cHead;
  93   chainTail = cTail;
  94   next = this;
  95   prev = this;
  96
  97   nextPolygon = NULL;
  98
  99   //compute bounding box
 100   directedLine* temp;
 101   minX = maxX = chainTail->head()[0];
 102   minY = maxY = chainTail->head()[1];
 103
 104   for(temp=chainHead; temp!=cTail; temp = temp->getNext())
 105     {
 106       if(temp->head()[0] < minX)
 107         minX = temp->head()[0];
 108       if(temp->head()[0] > maxX)
 109         maxX = temp->head()[0];
 110
 111       if(temp->head()[1] < minY)
 112         minY = temp->head()[1];
 113       if(temp->head()[1] > maxY)
 114         maxY = temp->head()[1];
 115     }
 116
 117   //check whether the chain is increasing or decreasing
 118   if(chainHead->compInY(chainTail) <0)
 119     isIncrease = 1;
 120   else
 121     isIncrease = 0;
 122
 123   //initilize currrent, this is used for accelerating search
 124   if(isIncrease)
 125     current = chainHead;
 126   else
 127     current = chainTail;
 128
 129   isKey = 0;
 130 }
 131
 132 //insert a new line between prev and this
 133 void monoChain::insert(monoChain* nc)
 134 {
 135   nc->next = this;
 136   nc->prev = prev;
 137   prev->next = nc;
 138   prev = nc;
 139 }
 140
 141 void monoChain::deleteLoop()
 142 {
 143   monoChain *temp, *tempNext;
 144   prev->next = NULL;
 145   for(temp=this; temp != NULL; temp = tempNext)
 146     {
 147       tempNext = temp->next;
 148       delete temp;
 149     }
 150 }
 151
 152 void monoChain::deleteLoopList()
 153 {
 154   monoChain *temp, *tempNext;
 155   for(temp=this; temp != NULL; temp = tempNext)
 156     {
 157       tempNext = temp->nextPolygon;
 158       temp->deleteLoop();
 159     }
 160 }
 161
 162 Int monoChain::toArraySingleLoop(monoChain** array, Int index)
 163 {
 164   monoChain *temp;
 165   array[index++] = this;
 166   for(temp = next; temp != this; temp = temp->next)
 167     {
 168       array[index++] = temp;
 169     }
 170   return index;
 171 }
 172
 173 monoChain** monoChain::toArrayAllLoops(Int& num_chains)
 174 {
 175   num_chains = numChainsAllLoops();
 176   monoChain **ret =  (monoChain**) malloc(sizeof(monoChain*) * num_chains);
 177   assert(ret);
 178   monoChain *temp;
 179   Int index = 0;
 180   for(temp = this; temp != NULL; temp=temp->nextPolygon){
 181     index = temp->toArraySingleLoop(ret, index);
 182   }
 183   return ret;
 184 }
 185
 186 Int monoChain::numChainsSingleLoop()
 187 {
 188   Int ret=0;
 189   monoChain* temp;
 190   if(next == this) return 1;
 191   ret = 1;
 192   for(temp=next; temp != this; temp = temp->next)
 193     ret++;
 194   return ret;
 195 }
 196
 197 Int monoChain::numChainsAllLoops()
 198 {
 199   Int ret=0;
 200   monoChain *temp;
 201   for(temp =this; temp != NULL; temp = temp->nextPolygon)
 202     ret += temp->numChainsSingleLoop();
 203   return ret;
 204 }
 205
 206 //update 'current'
 207 Real monoChain::chainIntersectHoriz(Real y)
 208 {
 209   directedLine* temp;
 210   if(isIncrease)
 211     {
 212       for(temp= current; temp != chainTail; temp = temp->getNext())
 213         {
 214           if(temp->head()[1] > y)
 215             break;
 216         }
 217       current = temp->getPrev();
 218     }
 219   else
 220     {
 221       for(temp = current; temp != chainHead; temp = temp->getPrev())
 222         {
 223           if(temp->head()[1] > y)
 224             break;
 225         }
 226       current = temp->getNext();
 227     }
 228   return intersectHoriz(current->head()[0], current->head()[1], current->tail()[0], current->tail()[1], y);
 229 }
 230
 231 monoChain* directedLineLoopToMonoChainLoop(directedLine* loop)
 232 {
 233   directedLine *temp;
 234   monoChain *ret=NULL;
 235
 236   //find the first cusp
 237   directedLine *prevCusp=NULL;
 238   directedLine *firstCusp;
 239
 240   if(isCusp(loop))
 241     prevCusp = loop;
 242   else
 243     {
 244       for(temp = loop->getNext(); temp != loop; temp = temp->getNext())
 245         if(isCusp(temp))
 246           break;
 247       prevCusp = temp;
 248     }
 249   firstCusp = prevCusp;
 250 //printf("first cusp is (%f,%f), (%f,%f), (%f,%f)\n", prevCusp->getPrev()->head()[0], prevCusp->getPrev()->head()[1], prevCusp->head()[0], prevCusp->head()[1], prevCusp->tail()[0], prevCusp->tail()[1]);
 251
 252   for(temp = prevCusp->getNext(); temp != loop; temp = temp->getNext())
 253     {
 254       if(isCusp(temp))
 255         {
 256 //printf("the cusp is (%f,%f), (%f,%f), (%f,%f)\n", temp->getPrev()->head()[0], temp->getPrev()->head()[1], temp->head()[0], temp->head()[1], temp->tail()[0], temp->tail()[1]);
 257           if(ret == NULL)
 258             {
 259               ret = new monoChain(prevCusp, temp);
 260             }
 261           else
 262             ret->insert(new monoChain(prevCusp, temp));
 263           prevCusp = temp;
 264         }
 265     }
 266   ret->insert(new monoChain(prevCusp, firstCusp));
 267
 268   return ret;
 269 }
 270
 271 monoChain* directedLineLoopListToMonoChainLoopList(directedLine* list)
 272 {
 273   directedLine* temp;
 274   monoChain* mc;
 275   monoChain* mcEnd;
 276   mc = directedLineLoopToMonoChainLoop(list);
 277   mcEnd = mc;
 278   for(temp = list->getNextPolygon(); temp != NULL; temp = temp->getNextPolygon())
 279     {
 280       monoChain *newLoop = directedLineLoopToMonoChainLoop(temp);
 281       mcEnd->setNextPolygon(newLoop);
 282       mcEnd = newLoop;
 283     }
 284   return mc;
 285 }
 286
 287 /*compare two edges of a polygon.
 288  *edge A < edge B if there is a horizontal line so that the intersection
 289  *with A is to the left of the intersection with B.
 290  *This function is used in sweepY for the dynamic search tree insertion to
 291  *order the edges.
 292  * Implementation: (x_1,y_1) and (x_2, y_2)
 293  */
 294 static Int compEdges(directedLine *e1, directedLine *e2)
 295 {
 296   Real* head1 = e1->head();
 297   Real* tail1 = e1->tail();
 298   Real* head2 = e2->head();
 299   Real* tail2 = e2->tail();
 300 /*
 301   Real h10 = head1[0];
 302   Real h11 = head1[1];
 303   Real t10 = tail1[0];
 304   Real t11 = tail1[1];
 305   Real h20 = head2[0];
 306   Real h21 = head2[1];
 307   Real t20 = tail2[0];
 308   Real t21 = tail2[1];
 309 */
 310   Real e1_Ymax, e1_Ymin, e2_Ymax, e2_Ymin;
 311 /*
 312   if(h11>t11) {
 313     e1_Ymax= h11;
 314     e1_Ymin= t11;
 315   }
 316   else{
 317     e1_Ymax = t11;
 318     e1_Ymin = h11;
 319   }
 320
 321   if(h21>t21) {
 322     e2_Ymax= h21;
 323     e2_Ymin= t21;
 324   }
 325   else{
 326     e2_Ymax = t21;
 327     e2_Ymin = h21;
 328   }
 329 */
 330
 331   if(head1[1]>tail1[1]) {
 332     e1_Ymax= head1[1];
 333     e1_Ymin= tail1[1];
 334   }
 335   else{
 336     e1_Ymax = tail1[1];
 337     e1_Ymin = head1[1];
 338   }
 339
 340   if(head2[1]>tail2[1]) {
 341     e2_Ymax= head2[1];
 342     e2_Ymin= tail2[1];
 343   }
 344   else{
 345     e2_Ymax = tail2[1];
 346     e2_Ymin = head2[1];
 347   }
 348
 349
 350   /*Real e1_Ymax = max(head1[1], tail1[1]);*/ /*max(e1->head()[1], e1->tail()[1]);*/
 351   /*Real e1_Ymin = min(head1[1], tail1[1]);*/ /*min(e1->head()[1], e1->tail()[1]);*/
 352   /*Real e2_Ymax = max(head2[1], tail2[1]);*/ /*max(e2->head()[1], e2->tail()[1]);*/
 353   /*Real e2_Ymin = min(head2[1], tail2[1]);*/ /*min(e2->head()[1], e2->tail()[1]);*/
 354
 355   Real Ymax = min(e1_Ymax, e2_Ymax);
 356   Real Ymin = max(e1_Ymin, e2_Ymin);
 357
 358   Real y = 0.5*(Ymax + Ymin);
 359
 360 /*  Real x1 = intersectHoriz(e1->head()[0], e1->head()[1], e1->tail()[0], e1->tail()[1], y);
 361   Real x2 = intersectHoriz(e2->head()[0], e2->head()[1], e2->tail()[0], e2->tail()[1], y);
 362 */
 363 /*
 364   Real x1 = intersectHoriz(h10, h11, t10, t11, y);
 365   Real x2 = intersectHoriz(h20, h21, t20, t21, y);
 366 */
 367   Real x1 = intersectHoriz(head1[0], head1[1], tail1[0], tail1[1], y);
 368   Real x2 = intersectHoriz(head2[0], head2[1], tail2[0], tail2[1], y);
 369
 370   if(x1<= x2) return -1;
 371   else return 1;
 372 }
 373
 374 Int  compChains(monoChain* mc1, monoChain* mc2)
 375 {
 376   Real y;
 377   assert(mc1->isKey || mc2->isKey);
 378   if(mc1->isKey)
 379     y = mc1->keyY;
 380   else
 381     y = mc2->keyY;
 382   directedLine *d1 = mc1->find(y);
 383   directedLine *d2 = mc2->find(y);
 384   mc2->find(y);
 385 //  Real x1 = mc1->chainIntersectHoriz(y);
 386 //  Real x2 = mc2->chainIntersectHoriz(y);
 387   return   compEdges(d1, d2);
 388 }
 389
 390 //this function modifies current for efficiency
 391 directedLine* monoChain::find(Real y)
 392 {
 393   directedLine *ret;
 394   directedLine *temp;
 395   assert(current->head()[1] <= y);
 396   if(isIncrease)
 397     {
 398       assert(chainTail->head()[1] >=y);
 399       for(temp=current; temp!=chainTail; temp = temp->getNext())
 400         {
 401           if(temp->head()[1] > y)
 402             break;
 403         }
 404       current = temp->getPrev();
 405       ret = current;
 406     }
 407   else
 408     {
 409       for(temp=current; temp != chainHead; temp = temp->getPrev())
 410         {
 411           if(temp->head()[1] > y)
 412             break;
 413         }
 414       current = temp->getNext();
 415       ret = temp;
 416     }
 417   return ret;
 418 }
 419
 420 void monoChain::printOneChain()
 421 {
 422   directedLine* temp;
 423   for(temp = chainHead; temp != chainTail; temp = temp->getNext())
 424     {
 425       printf("(%f,%f) ", temp->head()[0], temp->head()[1]);
 426     }
 427   printf("(%f,%f) \n", chainTail->head()[0], chainTail->head()[1]);
 428 }
 429
 430 void monoChain::printChainLoop()
 431 {
 432   monoChain* temp;
 433   this->printOneChain();
 434   for(temp = next; temp != this; temp = temp->next)
 435     {
 436       temp->printOneChain();
 437     }
 438   printf("\n");
 439 }
 440
 441 void monoChain::printAllLoops()
 442 {
 443   monoChain* temp;
 444   for(temp=this; temp != NULL; temp = temp->nextPolygon)
 445     temp->printChainLoop();
 446 }
 447
 448 //return 1 if error occures
 449 Int MC_sweepY(Int nVertices, monoChain** sortedVertices, sweepRange** ret_ranges)
 450 {
 451   Int i;
 452   Real keyY;
 453   Int errOccur=0;
 454 //printf("enter MC_sweepY\n");
 455 //printf("nVertices=%i\n", nVertices);
 456   /*for each vertex in the sorted list, update the binary search tree.
 457    *and store the range information for each vertex.
 458    */
 459   treeNode* searchTree = NULL;
 460 //printf("nVertices=%i\n", nVertices);
 461   for(i=0; i<nVertices; i++)
 462     {
 463       monoChain* vert = sortedVertices[i];
 464       keyY = vert->getHead()->head()[1]; //the sweep line
 465       directedLine *dline = vert->getHead();
 466       directedLine *dlinePrev = dline->getPrev();
 467       if(isBelow(dline, dline) && isBelow(dline, dlinePrev))
 468         {
 469 //printf("case 1\n");
 470           //this<v and prev < v
 471           //delete both edges
 472           vert->isKey = 1;
 473           vert->keyY = keyY;
 474           treeNode* thisNode = TreeNodeFind(searchTree, vert, (Int (*) (void *, void *))compChains);
 475           vert->isKey = 0;
 476
 477           vert->getPrev()->isKey = 1;
 478           vert->getPrev()->keyY = keyY;
 479           treeNode* prevNode = TreeNodeFind(searchTree, vert->getPrev(), (Int (*) (void *, void *))compChains);
 480           vert->getPrev()->isKey = 0;
 481
 482           if(cuspType(dline) == 1)//interior cusp
 483             {
 484
 485               treeNode* leftEdge = TreeNodePredecessor(prevNode);
 486               treeNode* rightEdge = TreeNodeSuccessor(thisNode);
 487               if(leftEdge == NULL ||  rightEdge == NULL)
 488                 {
 489                   errOccur = 1;
 490                   goto JUMP_HERE;
 491                 }
 492
 493               directedLine* leftEdgeDline = ((monoChain* ) leftEdge->key)->find(keyY);
 494
 495
 496
 497               directedLine* rightEdgeDline = ((monoChain* ) rightEdge->key)->find(keyY);
 498
 499               ret_ranges[i] = sweepRangeMake(leftEdgeDline, 1, rightEdgeDline, 1);
 500             }
 501           else /*exterior cusp*/
 502             {
 503               ret_ranges[i] = sweepRangeMake( dline, 1, dlinePrev, 1);
 504             }
 505
 506           searchTree = TreeNodeDeleteSingleNode(searchTree, thisNode);
 507           searchTree = TreeNodeDeleteSingleNode(searchTree, prevNode);
 508
 509         }
 510       else if(isAbove(dline, dline) && isAbove(dline, dlinePrev))
 511         {
 512 //printf("case 2\n");
 513           //insert both edges
 514           treeNode* thisNode = TreeNodeMake(vert);
 515           treeNode* prevNode = TreeNodeMake(vert->getPrev());
 516
 517           vert->isKey = 1;
 518           vert->keyY = keyY;
 519           searchTree = TreeNodeInsert(searchTree, thisNode, (Int (*) (void *, void *))compChains);
 520           vert->isKey = 0;
 521
 522           vert->getPrev()->isKey = 1;
 523           vert->getPrev()->keyY = keyY;
 524           searchTree = TreeNodeInsert(searchTree, prevNode, (Int (*) (void *, void *))compChains);
 525           vert->getPrev()->isKey = 0;
 526
 527           if(cuspType(dline) == 1) //interior cusp
 528             {
 529 //printf("cuspType is 1\n");
 530               treeNode* leftEdge = TreeNodePredecessor(thisNode);
 531               treeNode* rightEdge = TreeNodeSuccessor(prevNode);
 532               if(leftEdge == NULL || rightEdge == NULL)
 533                 {
 534                   errOccur = 1;
 535                   goto JUMP_HERE;
 536                 }
 537 //printf("leftEdge is %i, rightEdge is %i\n", leftEdge, rightEdge);
 538               directedLine* leftEdgeDline = ((monoChain*) leftEdge->key)->find(keyY);
 539               directedLine* rightEdgeDline = ((monoChain*) rightEdge->key)->find(keyY);
 540               ret_ranges[i] = sweepRangeMake( leftEdgeDline, 1, rightEdgeDline, 1);
 541             }
 542           else //exterior cusp
 543             {
 544 //printf("cuspType is not 1\n");
 545               ret_ranges[i] = sweepRangeMake(dlinePrev, 1, dline, 1);
 546             }
 547         }
 548       else
 549         {
 550 //printf("%i,%i\n", isAbove(dline, dline), isAbove(dline, dlinePrev));
 551           errOccur = 1;
 552           goto JUMP_HERE;
 553
 554           fprintf(stderr, "error in MC_sweepY\n");
 555           exit(1);
 556         }
 557     }
 558
 559  JUMP_HERE:
 560   //finally clean up space: delete  the search tree
 561   TreeNodeDeleteWholeTree(searchTree);
 562   return errOccur;
 563 }
 564
 565 void MC_findDiagonals(Int total_num_edges, monoChain** sortedVertices,
 566                    sweepRange** ranges, Int& num_diagonals,
 567                    directedLine** diagonal_vertices)
 568 {
 569   Int i,j,k;
 570   k=0;
 571   //reset 'current' of all the monoChains
 572   for(i=0; i<total_num_edges; i++)
 573     sortedVertices[i]->resetCurrent();
 574
 575   for(i=0; i<total_num_edges; i++)
 576     {
 577       directedLine* vert = sortedVertices[i]->getHead();
 578       directedLine* thisEdge = vert;
 579       directedLine* prevEdge = vert->getPrev();
 580       if(isBelow(vert, thisEdge) && isBelow(vert, prevEdge) && compEdges(prevEdge, thisEdge)<0)
 581         {
 582           //this is an upward interior cusp
 583           diagonal_vertices[k++] = vert;
 584
 585           directedLine* leftEdge = ranges[i]->left;
 586           directedLine* rightEdge = ranges[i]->right;
 587
 588           directedLine* leftVert = leftEdge;
 589           directedLine* rightVert = rightEdge->getNext();
 590           assert(leftVert->head()[1] >= vert->head()[1]);
 591           assert(rightVert->head()[1] >= vert->head()[1]);
 592           directedLine* minVert = (leftVert->head()[1] <= rightVert->head()[1])?leftVert:rightVert;
 593           Int found = 0;
 594           for(j=i+1; j<total_num_edges; j++)
 595             {
 596               if(sortedVertices[j]->getHead()->head()[1] > minVert->head()[1])
 597                 break;
 598
 599               if(sweepRangeEqual(ranges[i], ranges[j]))
 600                 {
 601                   found = 1;
 602                   break;
 603                 }
 604             }
 605
 606           if(found)
 607             diagonal_vertices[k++] = sortedVertices[j]->getHead();
 608           else
 609             diagonal_vertices[k++] = minVert;
 610         }
 611       else if(isAbove(vert, thisEdge) && isAbove(vert, prevEdge) && compEdges(prevEdge, thisEdge)>0)
 612         {
 613           //downward interior cusp
 614           diagonal_vertices[k++] = vert;
 615           directedLine* leftEdge = ranges[i]->left;
 616           directedLine* rightEdge = ranges[i]->right;
 617           directedLine* leftVert = leftEdge->getNext();
 618           directedLine* rightVert = rightEdge;
 619           assert(leftVert->head()[1] <= vert->head()[1]);
 620           assert(rightVert->head()[1] <= vert->head()[1]);
 621           directedLine* maxVert = (leftVert->head()[1] > rightVert->head()[1])? leftVert:rightVert;
 622           Int found=0;
 623           for(j=i-1; j>=0; j--)
 624             {
 625               if(sortedVertices[j]->getHead()->head()[1] < maxVert->head()[1])
 626                 break;
 627               if(sweepRangeEqual(ranges[i], ranges[j]))
 628                 {
 629                   found = 1;
 630                   break;
 631                 }
 632             }
 633           if(found)
 634             diagonal_vertices[k++] = sortedVertices[j]->getHead();
 635           else
 636             diagonal_vertices[k++] = maxVert;
 637         }
 638     }
 639   num_diagonals = k/2;
 640 }
 641
 642
 643
 644
 645 directedLine* MC_partitionY(directedLine *polygons, sampledLine **retSampledLines)
 646 {
 647 //printf("enter mc_partitionY\n");
 648   Int total_num_chains = 0;
 649   monoChain* loopList = directedLineLoopListToMonoChainLoopList(polygons);
 650   monoChain** array = loopList->toArrayAllLoops(total_num_chains);
 651
 652   if(total_num_chains<=2) //there is just one single monotone polygon
 653     {
 654       loopList->deleteLoopList();
 655       free(array);
 656       *retSampledLines = NULL;
 657       return polygons;
 658     }
 659
 660 //loopList->printAllLoops();
 661 //printf("total_num_chains=%i\n", total_num_chains);
 662   quicksort( (void**)array, 0, total_num_chains-1, (Int (*)(void*, void*))compChainHeadInY);
 663 //printf("after quicksort\n");
 664
 665   sweepRange** ranges = (sweepRange**)malloc(sizeof(sweepRange*) * (total_num_chains));
 666   assert(ranges);
 667
 668   if(MC_sweepY(total_num_chains, array, ranges))
 669     {
 670       loopList->deleteLoopList();
 671       free(array);
 672       *retSampledLines = NULL;
 673       return NULL;
 674     }
 675 //printf("after MC_sweepY\n");
 676
 677
 678   Int num_diagonals;
 679   /*number diagonals is < total_num_edges*total_num_edges*/
 680   directedLine** diagonal_vertices = (directedLine**) malloc(sizeof(directedLine*) * total_num_chains*2/*total_num_edges*/);
 681   assert(diagonal_vertices);
 682
 683 //printf("before call MC_findDiagonales\n");
 684
 685   MC_findDiagonals(total_num_chains, array, ranges, num_diagonals, diagonal_vertices);
 686 //printf("after call MC_findDia, num_diagnla=%i\n", num_diagonals);
 687
 688   directedLine* ret_polygons = polygons;
 689   sampledLine* newSampledLines = NULL;
 690   Int i,k;
 691
 692   num_diagonals=deleteRepeatDiagonals(num_diagonals, diagonal_vertices, diagonal_vertices);
 693
 694
 695
 696 //drawDiagonals(num_diagonals, diagonal_vertices);
 697 //printf("diagoanls are \n");
 698 //for(i=0; i<num_diagonals; i++)
 699 //  {
 700 //    printf("(%f,%f)\n", diagonal_vertices[2*i]->head()[0], diagonal_vertices[2*i]->head()[1]);
 701 //    printf("**(%f,%f)\n", diagonal_vertices[2*i+1]->head()[0], diagonal_vertices[2*i+1]->head()[1]);
 702 //  }
 703
 704   Int *removedDiagonals=(Int*)malloc(sizeof(Int) * num_diagonals);
 705   for(i=0; i<num_diagonals; i++)
 706     removedDiagonals[i] = 0;
 707 //  printf("first pass\n");
 708
 709
 710  for(i=0,k=0; i<num_diagonals; i++,k+=2)
 711     {
 712
 713
 714       directedLine* v1=diagonal_vertices[k];
 715       directedLine* v2=diagonal_vertices[k+1];
 716       directedLine* ret_p1;
 717       directedLine* ret_p2;
 718
 719       /*we ahve to determine whether v1 and v2 belong to the same polygon before
 720        *their structure are modified by connectDiagonal().
 721        */
 722 /*
 723       directedLine *root1 = v1->findRoot();
 724       directedLine *root2 = v2->findRoot();
 725       assert(root1);
 726       assert(root2);
 727 */
 728
 729 directedLine*  root1 = v1->rootLinkFindRoot();
 730 directedLine*  root2 = v2->rootLinkFindRoot();
 731
 732       if(root1 != root2)
 733         {
 734
 735           removedDiagonals[i] = 1;
 736           sampledLine* generatedLine;
 737
 738
 739
 740           v1->connectDiagonal(v1,v2, &ret_p1, &ret_p2, &generatedLine, ret_polygons);
 741
 742
 743
 744           newSampledLines = generatedLine->insert(newSampledLines);
 745 /*
 746           ret_polygons = ret_polygons->cutoffPolygon(root1);
 747
 748           ret_polygons = ret_polygons->cutoffPolygon(root2);
 749           ret_polygons = ret_p1->insertPolygon(ret_polygons);
 750 root1->rootLinkSet(ret_p1);
 751 root2->rootLinkSet(ret_p1);
 752 ret_p1->rootLinkSet(NULL);
 753 ret_p2->rootLinkSet(ret_p1);
 754 */
 755           ret_polygons = ret_polygons->cutoffPolygon(root2);
 756
 757
 758
 759 root2->rootLinkSet(root1);
 760 ret_p1->rootLinkSet(root1);
 761 ret_p2->rootLinkSet(root1);
 762
 763        /*now that we have connected the diagonal v1 and v2,
 764         *we have to check those unprocessed diagonals which
 765         *have v1 or v2 as an end point. Notice that the head of v1
 766         *has the same coodinates as the head of v2->prev, and the head of
 767         *v2 has the same coordinate as the head of v1->prev.
 768         *Suppose these is a diagonal (v1, x). If (v1,x) is still a valid
 769         *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
 770         *replaced by (v2->prev, x), that is, x is on the left of
 771         * v2->prev->prev->head, v2->prev->head, v2->prev->tail.
 772         */
 773         Int ii, kk;
 774         for(ii=0, kk=0; ii<num_diagonals; ii++, kk+=2)
 775           if( removedDiagonals[ii]==0)
 776             {
 777               directedLine* d1=diagonal_vertices[kk];
 778               directedLine* d2=diagonal_vertices[kk+1];
 779               /*check d1, and replace diagonal_vertices[kk] if necessary*/
 780               if(d1 == v1) {
 781                 /*check if d2 is to left of v1->prev->head:v1->head:v1->tail*/
 782                 if(! pointLeft2Lines(v1->getPrev()->head(),
 783                                      v1->head(), v1->tail(), d2->head()))
 784                   {
 785 /*
 786                     assert(pointLeft2Lines(v2->getPrev()->getPrev()->head(),
 787                                            v2->getPrev()->head(),
 788                                            v2->getPrev()->tail(), d2->head()));
 789 */
 790                     diagonal_vertices[kk] = v2->getPrev();
 791                   }
 792               }
 793               if(d1 == v2) {
 794                 /*check if d2 is to left of v2->prev->head:v2->head:v2->tail*/
 795                 if(! pointLeft2Lines(v2->getPrev()->head(),
 796                                      v2->head(), v2->tail(), d2->head()))
 797                   {
 798 /*
 799                     assert(pointLeft2Lines(v1->getPrev()->getPrev()->head(),
 800                                            v1->getPrev()->head(),
 801                                            v1->getPrev()->tail(), d2->head()));
 802 */
 803                     diagonal_vertices[kk] = v1->getPrev();
 804                   }
 805               }
 806               /*check d2 and replace diagonal_vertices[k+1] if necessary*/
 807               if(d2 == v1) {
 808                 /*check if d1 is to left of v1->prev->head:v1->head:v1->tail*/
 809                 if(! pointLeft2Lines(v1->getPrev()->head(),
 810                                      v1->head(), v1->tail(), d1->head()))
 811                   {
 812 /*                  assert(pointLeft2Lines(v2->getPrev()->getPrev()->head(),
 813                                            v2->getPrev()->head(),
 814                                            v2->getPrev()->tail(), d1->head()));
 815 */
 816                     diagonal_vertices[kk+1] = v2->getPrev();
 817                   }
 818               }
 819               if(d2 == v2) {
 820                 /*check if d1 is to left of v2->prev->head:v2->head:v2->tail*/
 821                 if(! pointLeft2Lines(v2->getPrev()->head(),
 822                                      v2->head(), v2->tail(), d1->head()))
 823                   {
 824 /*                  assert(pointLeft2Lines(v1->getPrev()->getPrev()->head(),
 825                                            v1->getPrev()->head(),
 826                                            v1->getPrev()->tail(), d1->head()));
 827 */
 828                     diagonal_vertices[kk+1] = v1->getPrev();
 829                   }
 830               }
 831             }
 832 }/*end if (root1 not equal to root 2)*/
 833 }
 834
 835   /*second pass,  now all diagoals should belong to the same polygon*/
 836 //printf("second pass: \n");
 837
 838 //  for(i=0; i<num_diagonals; i++)
 839 //    printf("%i ", removedDiagonals[i]);
 840
 841
 842   for(i=0,k=0; i<num_diagonals; i++, k += 2)
 843     if(removedDiagonals[i] == 0)
 844       {
 845
 846
 847         directedLine* v1=diagonal_vertices[k];
 848         directedLine* v2=diagonal_vertices[k+1];
 849
 850
 851
 852         directedLine* ret_p1;
 853         directedLine* ret_p2;
 854
 855         /*we ahve to determine whether v1 and v2 belong to the same polygon before
 856          *their structure are modified by connectDiagonal().
 857          */
 858         directedLine *root1 = v1->findRoot();
 859 /*
 860         directedLine *root2 = v2->findRoot();
 861
 862
 863
 864         assert(root1);
 865         assert(root2);
 866         assert(root1 == root2);
 867   */
 868         sampledLine* generatedLine;
 869
 870
 871
 872         v1->connectDiagonal(v1,v2, &ret_p1, &ret_p2, &generatedLine, ret_polygons);
 873         newSampledLines = generatedLine->insert(newSampledLines);
 874
 875         ret_polygons = ret_polygons->cutoffPolygon(root1);
 876
 877         ret_polygons = ret_p1->insertPolygon(ret_polygons);
 878
 879         ret_polygons = ret_p2->insertPolygon(ret_polygons);
 880
 881
 882
 883         for(Int j=i+1; j<num_diagonals; j++)
 884           {
 885             if(removedDiagonals[j] ==0)
 886               {
 887
 888                 directedLine* temp1=diagonal_vertices[2*j];
 889                 directedLine* temp2=diagonal_vertices[2*j+1];
 890                if(temp1==v1 || temp1==v2 || temp2==v1 || temp2==v2)
 891                 if(! temp1->samePolygon(temp1, temp2))
 892                   {
 893                     /*if temp1 and temp2 are in different polygons,
 894                      *then one of them must be v1 or v2.
 895                      */
 896
 897
 898
 899                     assert(temp1==v1 || temp1 == v2 || temp2==v1 || temp2 ==v2);
 900                     if(temp1==v1)
 901                       {
 902                         diagonal_vertices[2*j] = v2->getPrev();
 903                       }
 904                     if(temp2==v1)
 905                       {
 906                         diagonal_vertices[2*j+1] = v2->getPrev();
 907                       }
 908                     if(temp1==v2)
 909                       {
 910                         diagonal_vertices[2*j] = v1->getPrev();
 911                       }
 912                     if(temp2==v2)
 913                       {
 914                         diagonal_vertices[2*j+1] = v1->getPrev();
 915                       }
 916                   }
 917               }
 918           }
 919
 920       }
 921
 922
 923   //clean up
 924   loopList->deleteLoopList();
 925   free(array);
 926   free(ranges);
 927   free(diagonal_vertices);
 928   free(removedDiagonals);
 929
 930   *retSampledLines = newSampledLines;
 931   return ret_polygons;
 932 }
 933
 934