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