src/glu/sgi/libnurbs/nurbtess/monoTriangulation.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 "gluos.h"
  41 #include "glimports.h"
  42 #include "zlassert.h"
  43
  44 #include "monoTriangulation.h"
  45 #include "polyUtil.h" /*for area*/
  46 #include "partitionX.h"
  47 #include "monoPolyPart.h"
  48
  49
  50
  51 extern  directedLine*  polygonConvert(directedLine* polygon);
  52
  53 /*poly is NOT deleted
  54  */
  55 void monoTriangulationOpt(directedLine* poly, primStream* pStream)
  56 {
  57   Int n_cusps;
  58   Int n_edges = poly->numEdges();
  59   directedLine** cusps = (directedLine**) malloc(sizeof(directedLine*)*n_edges);
  60   assert(cusps);
  61   findInteriorCuspsX(poly, n_cusps, cusps);
  62   if(n_cusps ==0) //u monotine
  63     {
  64       monoTriangulationFun(poly, compV2InX, pStream);
  65     }
  66   else if(n_cusps == 1) // one interior cusp
  67     {
  68       directedLine* new_polygon = polygonConvert(cusps[0]);
  69       directedLine* other = findDiagonal_singleCuspX(new_polygon);
  70       //<other> should NOT be null unless there are self-intersecting
  71       //trim curves. In that case, we don't want to core dump, instead,
  72       //we triangulate anyway, and print out error message.
  73       if(other == NULL)
  74         {
  75           monoTriangulationFun(poly, compV2InX, pStream);
  76         }
  77       else
  78         {
  79           directedLine* ret_p1;
  80           directedLine* ret_p2;
  81
  82           new_polygon->connectDiagonal_2slines(new_polygon, other,
  83                                                &ret_p1,
  84                                                &ret_p2,
  85                                                new_polygon);
  86
  87           monoTriangulationFun(ret_p1, compV2InX, pStream);
  88           monoTriangulationFun(ret_p2, compV2InX, pStream);
  89
  90           ret_p1->deleteSinglePolygonWithSline();
  91           ret_p2->deleteSinglePolygonWithSline();
  92         }
  93     }
  94   else
  95     {
  96       //we need a general partitionX funtion (supposed to be in partitionX.C,
  97       //not implemented yet. XXX
  98       monoTriangulationFun(poly, compV2InY, pStream);
  99     }
 100
 101   free(cusps);
 102 }
 103
 104 void monoTriangulationRecOpt(Real* topVertex, Real* botVertex,
 105                              vertexArray* left_chain, Int left_current,
 106                              vertexArray* right_chain, Int right_current,
 107                              primStream* pStream)
 108 {
 109   Int i,j;
 110   Int n_left = left_chain->getNumElements();
 111   Int n_right = right_chain->getNumElements();
 112   if(left_current>= n_left-1 ||
 113      right_current>= n_right-1)
 114     {
 115       monoTriangulationRec(topVertex, botVertex, left_chain, left_current,
 116                            right_chain, right_current, pStream);
 117       return;
 118     }
 119   //now both left and right  have at least two vertices each.
 120   Real left_v = left_chain->getVertex(left_current)[1];
 121   Real right_v = right_chain->getVertex(right_current)[1];
 122
 123   if(left_v <= right_v) //first left vertex is below right
 124     {
 125       //find the last vertex of right which is above or equal to left
 126       for(j=right_current; j<=n_right-1; j++)
 127         {
 128           if(right_chain->getVertex(j)[1] < left_v)
 129             break;
 130         }
 131       monoTriangulationRecGen(topVertex, left_chain->getVertex(left_current),
 132                               left_chain, left_current, left_current,
 133                               right_chain, right_current, j-1,
 134                               pStream);
 135       monoTriangulationRecOpt(right_chain->getVertex(j-1),
 136                               botVertex,
 137                               left_chain, left_current,
 138                               right_chain, j,
 139                               pStream);
 140     }
 141   else //first right vertex is strictly below left
 142     {
 143       //find the last vertex of left which is strictly above right
 144       for(i=left_current; i<=n_left-1; i++)
 145         {
 146           if(left_chain->getVertex(i)[1] <= right_v)
 147             break;
 148         }
 149       monoTriangulationRecGen(topVertex, right_chain->getVertex(right_current),
 150                               left_chain, left_current, i-1,
 151                               right_chain, right_current, right_current,
 152                               pStream);
 153       monoTriangulationRecOpt(left_chain->getVertex(i-1),
 154                               botVertex,
 155                               left_chain, i,
 156                               right_chain, right_current,
 157                               pStream);
 158     }
 159 }
 160
 161
 162 void monoTriangulationRecGenTBOpt(Real* topVertex, Real* botVertex,
 163                           vertexArray* inc_chain, Int inc_current, Int inc_end,
 164                           vertexArray* dec_chain, Int dec_current, Int dec_end,
 165                           primStream* pStream)
 166 {
 167   pStream->triangle(topVertex, inc_chain->getVertex(inc_current), dec_chain->getVertex(dec_current));
 168
 169 /*printf("**(%f,%f)\n", inc_chain->getArray()[0][0],inc_chain->getArray()[0][1]);*/
 170   triangulateXYMonoTB(inc_end-inc_current+1, inc_chain->getArray()+inc_current,  dec_end-dec_current+1,  dec_chain->getArray()+dec_current, pStream);
 171
 172   pStream->triangle(botVertex, dec_chain->getVertex(dec_end), inc_chain->getVertex(inc_end));
 173 }
 174
 175
 176 /*n_left>=1
 177  *n_right>=1
 178  *the strip is going top to bottom. compared to the funtion
 179  * triangulateXYmono()
 180  */
 181 void triangulateXYMonoTB(Int n_left, Real** leftVerts,
 182                        Int n_right, Real** rightVerts,
 183                        primStream* pStream)
 184 {
 185
 186
 187   Int i,j,k,l;
 188   Real* topMostV;
 189
 190   assert(n_left>=1 && n_right>=1);
 191   if(leftVerts[0][1] >= rightVerts[0][1])
 192     {
 193       i=1;
 194       j=0;
 195       topMostV = leftVerts[0];
 196     }
 197   else
 198     {
 199       i=0;
 200       j=1;
 201       topMostV = rightVerts[0];
 202     }
 203
 204   while(1)
 205     {
 206       if(i >= n_left) /*case1: no more in left*/
 207         {
 208
 209           if(j<n_right-1) /*at least two vertices in right*/
 210             {
 211               pStream->begin();
 212               pStream->insert(topMostV);
 213               for(k=n_right-1; k>=j; k--)
 214                 pStream->insert(rightVerts[j]);
 215
 216               pStream->end(PRIMITIVE_STREAM_FAN);
 217
 218             }
 219
 220           break;
 221         }
 222       else if(j>= n_right) /*case2: no more in right*/
 223         {
 224
 225           if(i<n_left-1) /*at least two vertices in left*/
 226             {
 227               pStream->begin();
 228               pStream->insert(topMostV);
 229
 230               for(k=i; k<n_left; k++)
 231                 pStream->insert(leftVerts[k]);
 232
 233               pStream->end(PRIMITIVE_STREAM_FAN);
 234             }
 235
 236           break;
 237         }
 238       else /* case3: neither is empty, plus the topMostV, there is at least one triangle to output*/
 239         {
 240
 241           if(leftVerts[i][1] >=  rightVerts[j][1])
 242             {
 243               pStream->begin();
 244               pStream->insert(rightVerts[j]); /*the origin of this fan*/
 245
 246               pStream->insert(topMostV);
 247
 248               /*find the last k>=i such that
 249                *leftverts[k][1] >= rightverts[j][1]
 250                */
 251               k=i;
 252               while(k<n_left)
 253                 {
 254                   if(leftVerts[k][1] < rightVerts[j][1])
 255                     break;
 256                   k++;
 257                 }
 258               k--;
 259               for(l=i; l<=k; l++)
 260                 {
 261                   pStream->insert(leftVerts[l]);
 262                 }
 263
 264               pStream->end(PRIMITIVE_STREAM_FAN);
 265               //update i for next loop
 266               i = k+1;
 267               topMostV = leftVerts[k];
 268
 269             }
 270           else /*leftVerts[i][1] < rightVerts[j][1]*/
 271             {
 272               pStream->begin();
 273               pStream->insert(leftVerts[i]);/*the origion of this fan*/
 274
 275               /*find the last k>=j such that
 276                *rightverts[k][1] > leftverts[i][1]*/
 277               k=j;
 278               while(k< n_right)
 279                 {
 280                   if(rightVerts[k][1] <= leftVerts[i][1])
 281                     break;
 282                   k++;
 283                 }
 284               k--;
 285
 286               for(l=k; l>= j; l--)
 287                 pStream->insert(rightVerts[l]);
 288
 289               pStream->insert(topMostV);
 290               pStream->end(PRIMITIVE_STREAM_FAN);
 291               j=k+1;
 292               topMostV = rightVerts[j-1];
 293             }
 294         }
 295     }
 296 }
 297
 298 static int chainConvex(vertexArray* inc_chain, Int inc_current, Int inc_end)
 299 {
 300   Int i;
 301   //if there are no more than 2 vertices, return 1
 302   if(inc_current >= inc_end-1) return 1;
 303   for(i=inc_current; i<= inc_end-2; i++)
 304     {
 305       if(area(inc_chain->getVertex(i), inc_chain->getVertex(i+1), inc_chain->getVertex(i+2)) <0)
 306         return 0;
 307     }
 308   return 1;
 309 }
 310
 311 static int chainConcave(vertexArray* dec_chain, Int dec_current, Int dec_end)
 312 {
 313   Int i;
 314   //if there are no more than 2 vertices, return 1
 315   if(dec_current >= dec_end -1) return 1;
 316   for(i=dec_current; i<=dec_end-2; i++)
 317     {
 318       if(area(dec_chain->getVertex(i), dec_chain->getVertex(i+1), dec_chain->getVertex(i+2)) >0)
 319         return 0;
 320     }
 321   return 1;
 322 }
 323
 324 void monoTriangulationRecGenInU(Real* topVertex, Real* botVertex,
 325                                 vertexArray* inc_chain, Int inc_current, Int inc_end,
 326                                 vertexArray* dec_chain, Int dec_current, Int dec_end,
 327                                 primStream* pStream)
 328 {
 329
 330 }
 331
 332 void  monoTriangulationRecGenOpt(Real* topVertex, Real* botVertex,
 333                                  vertexArray* inc_chain, Int inc_current, Int inc_end,
 334                                  vertexArray* dec_chain, Int dec_current, Int dec_end,
 335                                  primStream* pStream)
 336 {
 337   Int i;
 338   //copy this to a polygon: directedLine Lioop
 339   sampledLine* sline;
 340   directedLine* dline;
 341   directedLine* poly;
 342
 343   if(inc_current <= inc_end) //at least one vertex in inc_chain
 344     {
 345       sline = new sampledLine(topVertex, inc_chain->getVertex(inc_current));
 346       poly = new directedLine(INCREASING, sline);
 347       for(i=inc_current; i<=inc_end-1; i++)
 348         {
 349           sline = new sampledLine(inc_chain->getVertex(i), inc_chain->getVertex(i+1));
 350           dline = new directedLine(INCREASING, sline);
 351           poly->insert(dline);
 352         }
 353       sline = new sampledLine(inc_chain->getVertex(inc_end), botVertex);
 354       dline = new directedLine(INCREASING, sline);
 355       poly->insert(dline);
 356     }
 357   else //inc_chian is empty
 358     {
 359       sline = new sampledLine(topVertex, botVertex);
 360       dline = new directedLine(INCREASING, sline);
 361       poly = dline;
 362     }
 363
 364   assert(poly != NULL);
 365
 366   if(dec_current <= dec_end) //at least on vertex in dec_Chain
 367     {
 368       sline = new sampledLine(botVertex, dec_chain->getVertex(dec_end));
 369       dline = new directedLine(INCREASING, sline);
 370       poly->insert(dline);
 371       for(i=dec_end; i>dec_current; i--)
 372         {
 373           sline = new sampledLine(dec_chain->getVertex(i), dec_chain->getVertex(i-1));
 374           dline = new directedLine(INCREASING, sline);
 375           poly->insert(dline);
 376         }
 377       sline = new sampledLine(dec_chain->getVertex(dec_current), topVertex);
 378       dline = new directedLine(INCREASING, sline);
 379       poly->insert(dline);
 380     }
 381   else //dec_chain  is empty
 382     {
 383       sline = new sampledLine(botVertex, topVertex);
 384       dline = new directedLine(INCREASING, sline);
 385       poly->insert(dline);
 386     }
 387
 388   {
 389     Int n_cusps;
 390     Int n_edges = poly->numEdges();
 391     directedLine** cusps = (directedLine**) malloc(sizeof(directedLine*)*n_edges);
 392     assert(cusps);
 393     findInteriorCuspsX(poly, n_cusps, cusps);
 394
 395     if(n_cusps ==0) //u monotine
 396       {
 397         monoTriangulationFun(poly, compV2InX, pStream);
 398       }
 399     else if(n_cusps == 1) // one interior cusp
 400       {
 401         directedLine* new_polygon = polygonConvert(cusps[0]);
 402         directedLine* other = findDiagonal_singleCuspX(new_polygon);
 403           //<other> should NOT be null unless there are self-intersecting
 404           //trim curves. In that case, we don't want to core dump, instead,
 405           //we triangulate anyway, and print out error message.
 406           if(other == NULL)
 407             {
 408               monoTriangulationFun(poly, compV2InX, pStream);
 409             }
 410           else
 411             {
 412               directedLine* ret_p1;
 413               directedLine* ret_p2;
 414
 415               new_polygon->connectDiagonal_2slines(new_polygon, other,
 416                                                    &ret_p1,
 417                                                    &ret_p2,
 418                                                    new_polygon);
 419
 420               monoTriangulationFun(ret_p1, compV2InX, pStream);
 421               monoTriangulationFun(ret_p2, compV2InX, pStream);
 422
 423               ret_p1->deleteSinglePolygonWithSline();
 424               ret_p2->deleteSinglePolygonWithSline();
 425             }
 426       }
 427     else
 428       {
 429         //we need a general partitionX funtion (supposed to be in partitionX.C,
 430         //not implemented yet. XXX
 431         //monoTriangulationFun(poly, compV2InY, pStream);
 432
 433         directedLine* new_polygon = polygonConvert(poly);
 434         directedLine* list = monoPolyPart(new_polygon);
 435         for(directedLine* temp = list; temp != NULL; temp = temp->getNextPolygon())
 436           {
 437             monoTriangulationFun(temp, compV2InX, pStream);
 438           }
 439         //clean up
 440         list->deletePolygonListWithSline();
 441
 442       }
 443
 444     free(cusps);
 445     /*
 446       if(numInteriorCuspsX(poly) == 0) //is u monotone
 447         monoTriangulationFun(poly, compV2InX, pStream);
 448       else //it is not u motone
 449         monoTriangulationFun(poly, compV2InY, pStream);
 450         */
 451     //clean up space
 452     poly->deleteSinglePolygonWithSline();
 453     return;
 454   }
 455
 456   //apparently the following code is not reachable,
 457   //it is for test purpose
 458   if(inc_current > inc_end || dec_current>dec_end)
 459     {
 460     monoTriangulationRecGen(topVertex, botVertex, inc_chain, inc_current, inc_end,
 461                             dec_chain, dec_current, dec_end,
 462                             pStream);
 463     return;
 464   }
 465
 466
 467   if(
 468      area(dec_chain->getVertex(dec_current),
 469           topVertex,
 470           inc_chain->getVertex(inc_current)) >=0
 471      && chainConvex(inc_chain, inc_current, inc_end)
 472      && chainConcave(dec_chain, dec_current, dec_end)
 473      && area(inc_chain->getVertex(inc_end), botVertex, dec_chain->getVertex(dec_end)) >=0
 474      )
 475     {
 476       monoTriangulationRecFunGen(topVertex, botVertex,
 477                                  inc_chain, inc_current, inc_end,
 478                                  dec_chain, dec_current, dec_end,
 479                                  compV2InX, pStream);
 480     }
 481   else
 482     {
 483       monoTriangulationRecGen(topVertex, botVertex, inc_chain, inc_current, inc_end,
 484                               dec_chain, dec_current, dec_end,
 485                               pStream);
 486     }
 487 }
 488
 489 /*if inc_current>inc_end, then inc_chain has no points to be considered
 490  *same for dec_chain
 491  */
 492 void monoTriangulationRecGen(Real* topVertex, Real* botVertex,
 493                           vertexArray* inc_chain, Int inc_current, Int inc_end,
 494                           vertexArray* dec_chain, Int dec_current, Int dec_end,
 495                           primStream* pStream)
 496 {
 497   Real** inc_array ;
 498   Real** dec_array ;
 499   Int i;
 500
 501   if(inc_current > inc_end && dec_current>dec_end)
 502     return;
 503   else if(inc_current>inc_end) /*no more vertices on inc_chain*/
 504     {
 505       dec_array = dec_chain->getArray();
 506       reflexChain rChain(100,0);
 507       /*put the top vertex into the reflex chain*/
 508       rChain.processNewVertex(topVertex, pStream);
 509       /*process all the vertices on the dec_chain*/
 510       for(i=dec_current; i<=dec_end; i++){
 511         rChain.processNewVertex(dec_array[i], pStream);
 512       }
 513       /*process the bottom vertex*/
 514       rChain.processNewVertex(botVertex, pStream);
 515     }
 516   else if(dec_current> dec_end) /*no more vertices on dec_chain*/
 517     {
 518       inc_array = inc_chain->getArray();
 519
 520       reflexChain rChain(100,1);
 521       /*put the top vertex into the reflex chain*/
 522       rChain.processNewVertex(topVertex, pStream);
 523       /*process all the vertices on the inc_chain*/
 524       for(i=inc_current; i<=inc_end; i++){
 525         rChain.processNewVertex(inc_array[i], pStream);
 526       }
 527       /*process the bottom vertex*/
 528       rChain.processNewVertex(botVertex, pStream);
 529     }
 530   else /*neither chain is empty*/
 531     {
 532       inc_array = inc_chain -> getArray();
 533       dec_array = dec_chain -> getArray();
 534
 535       /*if top of inc_chain is 'lower' than top of dec_chain, process all the
 536        *vertices on the dec_chain which are higher than top of inc_chain
 537        */
 538       if(compV2InY(inc_array[inc_current], dec_array[dec_current]) <= 0)
 539         {
 540
 541           reflexChain rChain(100, 0);
 542           rChain.processNewVertex(topVertex, pStream);
 543           for(i=dec_current; i<=dec_end; i++)
 544             {
 545               if(compV2InY(inc_array[inc_current], dec_array[i]) <= 0)
 546                 rChain.processNewVertex(dec_array[i], pStream);
 547               else
 548                 break;
 549             }
 550           rChain.outputFan(inc_array[inc_current], pStream);
 551           monoTriangulationRecGen(dec_array[i-1], botVertex,
 552                                inc_chain, inc_current, inc_end,
 553                                dec_chain, i, dec_end,
 554                                pStream);
 555         }
 556       else /*compV2InY(inc_array[inc_current], dec_array[dec_current]) > 0*/
 557         {
 558
 559           reflexChain rChain(100, 1);
 560           rChain.processNewVertex(topVertex, pStream);
 561           for(i=inc_current; i<=inc_end; i++)
 562             {
 563               if(compV2InY(inc_array[i], dec_array[dec_current]) >0)
 564                 rChain.processNewVertex(inc_array[i], pStream);
 565               else
 566                 break;
 567             }
 568           rChain.outputFan(dec_array[dec_current], pStream);
 569           monoTriangulationRecGen(inc_array[i-1], botVertex,
 570                                inc_chain, i, inc_end,
 571                                dec_chain, dec_current,dec_end,
 572                                pStream);
 573         }
 574     }/*end case neither is empty*/
 575 }
 576
 577 void monoTriangulationFun(directedLine* monoPolygon, Int (*compFun)(Real*, Real*), primStream* pStream)
 578 {
 579   Int i;
 580   /*find the top vertex, bottom vertex, inccreasing chain, and decreasing chain,
 581    *then call monoTriangulationRec
 582    */
 583   directedLine* tempV;
 584   directedLine* topV;
 585   directedLine* botV;
 586   topV = botV = monoPolygon;
 587   for(tempV = monoPolygon->getNext(); tempV != monoPolygon; tempV = tempV->getNext())
 588     {
 589       if(compFun(topV->head(), tempV->head())<0) {
 590         topV = tempV;
 591       }
 592       if(compFun(botV->head(), tempV->head())>0) {
 593         botV = tempV;
 594       }
 595     }
 596
 597   /*creat increase and decrease chains*/
 598   vertexArray inc_chain(20); /*this is a dynamic array*/
 599   for(i=1; i<=topV->get_npoints()-2; i++) { /*the first vertex is the top vertex which doesn't belong to inc_chain*/
 600     inc_chain.appendVertex(topV->getVertex(i));
 601   }
 602   for(tempV = topV->getNext(); tempV != botV; tempV = tempV->getNext())
 603     {
 604       for(i=0; i<=tempV->get_npoints()-2; i++){
 605         inc_chain.appendVertex(tempV->getVertex(i));
 606       }
 607     }
 608
 609   vertexArray dec_chain(20);
 610   for(tempV = topV->getPrev(); tempV != botV; tempV = tempV->getPrev())
 611     {
 612       for(i=tempV->get_npoints()-2; i>=0; i--){
 613         dec_chain.appendVertex(tempV->getVertex(i));
 614       }
 615     }
 616   for(i=botV->get_npoints()-2; i>=1; i--){
 617     dec_chain.appendVertex(tempV->getVertex(i));
 618   }
 619
 620   if (!(0 == inc_chain.getNumElements() && 0 == dec_chain.getNumElements())) {
 621      monoTriangulationRecFun(topV->head(), botV->head(), &inc_chain, 0,
 622                              &dec_chain, 0, compFun, pStream);
 623   }
 624 }
 625
 626 void monoTriangulation(directedLine* monoPolygon, primStream* pStream)
 627 {
 628   Int i;
 629   /*find the top vertex, bottom vertex, inccreasing chain, and decreasing chain,
 630    *then call monoTriangulationRec
 631    */
 632   directedLine* tempV;
 633   directedLine* topV;
 634   directedLine* botV;
 635   topV = botV = monoPolygon;
 636   for(tempV = monoPolygon->getNext(); tempV != monoPolygon; tempV = tempV->getNext())
 637     {
 638       if(compV2InY(topV->head(), tempV->head())<0) {
 639         topV = tempV;
 640       }
 641       if(compV2InY(botV->head(), tempV->head())>0) {
 642         botV = tempV;
 643       }
 644     }
 645   /*creat increase and decrease chains*/
 646   vertexArray inc_chain(20); /*this is a dynamic array*/
 647   for(i=1; i<=topV->get_npoints()-2; i++) { /*the first vertex is the top vertex which doesn't belong to inc_chain*/
 648     inc_chain.appendVertex(topV->getVertex(i));
 649   }
 650   for(tempV = topV->getNext(); tempV != botV; tempV = tempV->getNext())
 651     {
 652       for(i=0; i<=tempV->get_npoints()-2; i++){
 653         inc_chain.appendVertex(tempV->getVertex(i));
 654       }
 655     }
 656
 657   vertexArray dec_chain(20);
 658   for(tempV = topV->getPrev(); tempV != botV; tempV = tempV->getPrev())
 659     {
 660       for(i=tempV->get_npoints()-2; i>=0; i--){
 661         dec_chain.appendVertex(tempV->getVertex(i));
 662       }
 663     }
 664   for(i=botV->get_npoints()-2; i>=1; i--){
 665     dec_chain.appendVertex(tempV->getVertex(i));
 666   }
 667
 668   monoTriangulationRec(topV->head(), botV->head(), &inc_chain, 0, &dec_chain, 0, pStream);
 669
 670 }
 671
 672 /*the chain could be increasing or decreasing, although we use the
 673  * name inc_chain.
 674  *the argument  is_increase_chain indicates whether this chain
 675  *is increasing (left chain in V-monotone case) or decreaing (right chain
 676  *in V-monotone case).
 677  */
 678 void monoTriangulation2(Real* topVertex, Real* botVertex,
 679                         vertexArray* inc_chain, Int inc_smallIndex,
 680                         Int inc_largeIndex,
 681                         Int is_increase_chain,
 682                         primStream* pStream)
 683 {
 684   assert( inc_chain != NULL);
 685   Real** inc_array ;
 686
 687   if(inc_smallIndex > inc_largeIndex)
 688     return; //no triangles
 689   if(inc_smallIndex == inc_largeIndex)
 690     {
 691       if(is_increase_chain)
 692         pStream->triangle(inc_chain->getVertex(inc_smallIndex), botVertex, topVertex);
 693       else
 694         pStream->triangle(inc_chain->getVertex(inc_smallIndex), topVertex, botVertex);
 695       return;
 696     }
 697   Int i;
 698
 699   if(is_increase_chain && botVertex[1] == inc_chain->getVertex(inc_largeIndex)[1])
 700     {
 701       pStream->triangle(botVertex, inc_chain->getVertex(inc_largeIndex-1),
 702                         inc_chain->getVertex(inc_largeIndex));
 703       monoTriangulation2(topVertex, botVertex, inc_chain, inc_smallIndex,
 704                          inc_largeIndex-1,
 705                          is_increase_chain,
 706                          pStream);
 707       return;
 708     }
 709   else if( (!is_increase_chain) && topVertex[1] == inc_chain->getVertex(inc_smallIndex)[1])
 710     {
 711       pStream->triangle(topVertex, inc_chain->getVertex(inc_smallIndex+1),
 712                         inc_chain->getVertex(inc_smallIndex));
 713       monoTriangulation2(topVertex, botVertex, inc_chain, inc_smallIndex+1,
 714                          inc_largeIndex, is_increase_chain, pStream);
 715       return ;
 716     }
 717
 718   inc_array = inc_chain->getArray();
 719
 720   reflexChain rChain(20,is_increase_chain); /*1 means the chain is increasing*/
 721
 722   rChain.processNewVertex(topVertex, pStream);
 723
 724   for(i=inc_smallIndex; i<=inc_largeIndex; i++){
 725     rChain.processNewVertex(inc_array[i], pStream);
 726   }
 727   rChain.processNewVertex(botVertex, pStream);
 728
 729 }
 730
 731 /*if compFun == compV2InY, top to bottom: V-monotone
 732  *if compFun == compV2InX, right to left: U-monotone
 733  */
 734 void monoTriangulationRecFunGen(Real* topVertex, Real* botVertex,
 735                           vertexArray* inc_chain, Int inc_current, Int inc_end,
 736                           vertexArray* dec_chain, Int dec_current, Int dec_end,
 737                           Int  (*compFun)(Real*, Real*),
 738                           primStream* pStream)
 739 {
 740   assert( inc_chain != NULL && dec_chain != NULL);
 741   assert( ! (inc_current> inc_end &&
 742              dec_current> dec_end));
 743   /*
 744   Int inc_nVertices;
 745   Int dec_nVertices;
 746   */
 747   Real** inc_array ;
 748   Real** dec_array ;
 749   Int i;
 750   assert( ! ( (inc_chain==NULL) && (dec_chain==NULL)));
 751
 752   if(inc_current> inc_end) /*no more vertices on inc_chain*/
 753     {
 754
 755       dec_array = dec_chain->getArray();
 756       reflexChain rChain(20,0);
 757       /*put the top vertex into the reflex chain*/
 758       rChain.processNewVertex(topVertex, pStream);
 759       /*process all the vertices on the dec_chain*/
 760       for(i=dec_current; i<=dec_end; i++){
 761         rChain.processNewVertex(dec_array[i], pStream);
 762       }
 763       /*process the bottom vertex*/
 764       rChain.processNewVertex(botVertex, pStream);
 765
 766     }
 767   else if(dec_current> dec_end) /*no more vertices on dec_chain*/
 768     {
 769       inc_array = inc_chain->getArray();
 770       reflexChain rChain(20,1);
 771       /*put the top vertex into the reflex chain*/
 772       rChain.processNewVertex(topVertex, pStream);
 773       /*process all the vertices on the inc_chain*/
 774       for(i=inc_current; i<=inc_end; i++){
 775         rChain.processNewVertex(inc_array[i], pStream);
 776       }
 777       /*process the bottom vertex*/
 778       rChain.processNewVertex(botVertex, pStream);
 779     }
 780   else /*neither chain is empty*/
 781     {
 782       inc_array = inc_chain -> getArray();
 783       dec_array = dec_chain -> getArray();
 784
 785       /*if top of inc_chain is 'lower' than top of dec_chain, process all the
 786        *vertices on the dec_chain which are higher than top of inc_chain
 787        */
 788       if(compFun(inc_array[inc_current], dec_array[dec_current]) <= 0)
 789         {
 790
 791           reflexChain rChain(20, 0);
 792           rChain.processNewVertex(topVertex, pStream);
 793           for(i=dec_current; i<=dec_end; i++)
 794             {
 795               if(compFun(inc_array[inc_current], dec_array[i]) <= 0)
 796                 rChain.processNewVertex(dec_array[i], pStream);
 797               else
 798                 break;
 799             }
 800           rChain.outputFan(inc_array[inc_current], pStream);
 801           monoTriangulationRecFunGen(dec_array[i-1], botVertex,
 802                                inc_chain, inc_current, inc_end,
 803                                dec_chain, i, dec_end,
 804                                compFun,
 805                                pStream);
 806         }
 807       else /*compFun(inc_array[inc_current], dec_array[dec_current]) > 0*/
 808         {
 809
 810           reflexChain rChain(20, 1);
 811           rChain.processNewVertex(topVertex, pStream);
 812           for(i=inc_current; i<=inc_end; i++)
 813             {
 814               if(compFun(inc_array[i], dec_array[dec_current]) >0)
 815                 rChain.processNewVertex(inc_array[i], pStream);
 816               else
 817                 break;
 818             }
 819           rChain.outputFan(dec_array[dec_current], pStream);
 820           monoTriangulationRecFunGen(inc_array[i-1], botVertex,
 821                                inc_chain, i,inc_end,
 822                                dec_chain, dec_current,dec_end,
 823                                compFun,
 824                                pStream);
 825         }
 826     }/*end case neither is empty*/
 827 }
 828
 829 /*if compFun == compV2InY, top to bottom: V-monotone
 830  *if compFun == compV2InX, right to left: U-monotone
 831  */
 832 void monoTriangulationRecFun(Real* topVertex, Real* botVertex,
 833                           vertexArray* inc_chain, Int inc_current,
 834                           vertexArray* dec_chain, Int dec_current,
 835                           Int  (*compFun)(Real*, Real*),
 836                           primStream* pStream)
 837 {
 838   assert( inc_chain != NULL && dec_chain != NULL);
 839   assert( ! (inc_current>=inc_chain->getNumElements() &&
 840              dec_current>=dec_chain->getNumElements()));
 841   Int inc_nVertices;
 842   Int dec_nVertices;
 843   Real** inc_array ;
 844   Real** dec_array ;
 845   Int i;
 846   assert( ! ( (inc_chain==NULL) && (dec_chain==NULL)));
 847
 848   if(inc_current>=inc_chain->getNumElements()) /*no more vertices on inc_chain*/
 849     {
 850
 851       dec_array = dec_chain->getArray();
 852       dec_nVertices = dec_chain->getNumElements();
 853       reflexChain rChain(20,0);
 854       /*put the top vertex into the reflex chain*/
 855       rChain.processNewVertex(topVertex, pStream);
 856       /*process all the vertices on the dec_chain*/
 857       for(i=dec_current; i<dec_nVertices; i++){
 858         rChain.processNewVertex(dec_array[i], pStream);
 859       }
 860       /*process the bottom vertex*/
 861       rChain.processNewVertex(botVertex, pStream);
 862
 863     }
 864   else if(dec_current>= dec_chain->getNumElements()) /*no more vertices on dec_chain*/
 865     {
 866       inc_array = inc_chain->getArray();
 867       inc_nVertices= inc_chain->getNumElements();
 868       reflexChain rChain(20,1);
 869       /*put the top vertex into the reflex chain*/
 870       rChain.processNewVertex(topVertex, pStream);
 871       /*process all the vertices on the inc_chain*/
 872       for(i=inc_current; i<inc_nVertices; i++){
 873         rChain.processNewVertex(inc_array[i], pStream);
 874       }
 875       /*process the bottom vertex*/
 876       rChain.processNewVertex(botVertex, pStream);
 877     }
 878   else /*neither chain is empty*/
 879     {
 880       inc_array = inc_chain -> getArray();
 881       dec_array = dec_chain -> getArray();
 882       inc_nVertices= inc_chain->getNumElements();
 883       dec_nVertices= dec_chain->getNumElements();
 884       /*if top of inc_chain is 'lower' than top of dec_chain, process all the
 885        *vertices on the dec_chain which are higher than top of inc_chain
 886        */
 887       if(compFun(inc_array[inc_current], dec_array[dec_current]) <= 0)
 888         {
 889
 890           reflexChain rChain(20, 0);
 891           rChain.processNewVertex(topVertex, pStream);
 892           for(i=dec_current; i<dec_nVertices; i++)
 893             {
 894               if(compFun(inc_array[inc_current], dec_array[i]) <= 0)
 895                 rChain.processNewVertex(dec_array[i], pStream);
 896               else
 897                 break;
 898             }
 899           rChain.outputFan(inc_array[inc_current], pStream);
 900           monoTriangulationRecFun(dec_array[i-1], botVertex,
 901                                inc_chain, inc_current,
 902                                dec_chain, i,
 903                                compFun,
 904                                pStream);
 905         }
 906       else /*compFun(inc_array[inc_current], dec_array[dec_current]) > 0*/
 907         {
 908
 909           reflexChain rChain(20, 1);
 910           rChain.processNewVertex(topVertex, pStream);
 911           for(i=inc_current; i<inc_nVertices; i++)
 912             {
 913               if(compFun(inc_array[i], dec_array[dec_current]) >0)
 914                 rChain.processNewVertex(inc_array[i], pStream);
 915               else
 916                 break;
 917             }
 918           rChain.outputFan(dec_array[dec_current], pStream);
 919           monoTriangulationRecFun(inc_array[i-1], botVertex,
 920                                inc_chain, i,
 921                                dec_chain, dec_current,
 922                                compFun,
 923                                pStream);
 924         }
 925     }/*end case neither is empty*/
 926 }
 927
 928
 929 void monoTriangulationRec(Real* topVertex, Real* botVertex,
 930                           vertexArray* inc_chain, Int inc_current,
 931                           vertexArray* dec_chain, Int dec_current,
 932                           primStream* pStream)
 933 {
 934   assert( inc_chain != NULL && dec_chain != NULL);
 935   assert( ! (inc_current>=inc_chain->getNumElements() &&
 936              dec_current>=dec_chain->getNumElements()));
 937   Int inc_nVertices;
 938   Int dec_nVertices;
 939   Real** inc_array ;
 940   Real** dec_array ;
 941   Int i;
 942   assert( ! ( (inc_chain==NULL) && (dec_chain==NULL)));
 943
 944   if(inc_current>=inc_chain->getNumElements()) /*no more vertices on inc_chain*/
 945     {
 946
 947       dec_array = dec_chain->getArray();
 948       dec_nVertices = dec_chain->getNumElements();
 949       reflexChain rChain(20,0);
 950       /*put the top vertex into the reflex chain*/
 951       rChain.processNewVertex(topVertex, pStream);
 952       /*process all the vertices on the dec_chain*/
 953       for(i=dec_current; i<dec_nVertices; i++){
 954         rChain.processNewVertex(dec_array[i], pStream);
 955       }
 956       /*process the bottom vertex*/
 957       rChain.processNewVertex(botVertex, pStream);
 958
 959     }
 960   else if(dec_current>= dec_chain->getNumElements()) /*no more vertices on dec_chain*/
 961     {
 962       inc_array = inc_chain->getArray();
 963       inc_nVertices= inc_chain->getNumElements();
 964       reflexChain rChain(20,1);
 965       /*put the top vertex into the reflex chain*/
 966       rChain.processNewVertex(topVertex, pStream);
 967       /*process all the vertices on the inc_chain*/
 968       for(i=inc_current; i<inc_nVertices; i++){
 969         rChain.processNewVertex(inc_array[i], pStream);
 970       }
 971       /*process the bottom vertex*/
 972       rChain.processNewVertex(botVertex, pStream);
 973     }
 974   else /*neither chain is empty*/
 975     {
 976       inc_array = inc_chain -> getArray();
 977       dec_array = dec_chain -> getArray();
 978       inc_nVertices= inc_chain->getNumElements();
 979       dec_nVertices= dec_chain->getNumElements();
 980       /*if top of inc_chain is 'lower' than top of dec_chain, process all the
 981        *vertices on the dec_chain which are higher than top of inc_chain
 982        */
 983       if(compV2InY(inc_array[inc_current], dec_array[dec_current]) <= 0)
 984         {
 985
 986           reflexChain rChain(20, 0);
 987           rChain.processNewVertex(topVertex, pStream);
 988           for(i=dec_current; i<dec_nVertices; i++)
 989             {
 990               if(compV2InY(inc_array[inc_current], dec_array[i]) <= 0)
 991                 rChain.processNewVertex(dec_array[i], pStream);
 992               else
 993                 break;
 994             }
 995           rChain.outputFan(inc_array[inc_current], pStream);
 996           monoTriangulationRec(dec_array[i-1], botVertex,
 997                                inc_chain, inc_current,
 998                                dec_chain, i,
 999                                pStream);
1000         }
1001       else /*compV2InY(inc_array[inc_current], dec_array[dec_current]) > 0*/
1002         {
1003
1004           reflexChain rChain(20, 1);
1005           rChain.processNewVertex(topVertex, pStream);
1006           for(i=inc_current; i<inc_nVertices; i++)
1007             {
1008               if(compV2InY(inc_array[i], dec_array[dec_current]) >0)
1009                 rChain.processNewVertex(inc_array[i], pStream);
1010               else
1011                 break;
1012             }
1013           rChain.outputFan(dec_array[dec_current], pStream);
1014           monoTriangulationRec(inc_array[i-1], botVertex,
1015                                inc_chain, i,
1016                                dec_chain, dec_current,
1017                                pStream);
1018         }
1019     }/*end case neither is empty*/
1020 }
1021
1022
1023
1024 /* the name here assumes that the polygon is Y-monotone, but
1025  *this function also works for X-monotone polygons.
1026  * a monotne polygon consists of two extrem verteices: topVertex and botVertex, and
1027  *two monotone chains: inc_chain, and dec_chain. The edges of the increasing chain (inc_chain)
1028  *is ordered by following pointer: next, while  the edges of the decreasing chain (dec_chain)
1029  *is ordered by following pointer: prev
1030  * inc_index index the vertex which is the toppest of the inc_chain which we are handling currently.
1031  * dec_index index the vertex which is the toppest of the dec_chain which we are handling currently.
1032  */
1033 void monoTriangulationRec(directedLine* inc_chain, Int inc_index,
1034                           directedLine* dec_chain, Int dec_index,
1035                           directedLine* topVertex, Int top_index,
1036                           directedLine* botVertex,
1037                           primStream* pStream)
1038 {
1039   Int i;
1040   directedLine *temp, *oldtemp = NULL;
1041   Int tempIndex, oldtempIndex = 0;
1042
1043   assert(inc_chain != NULL && dec_chain != NULL);
1044
1045   if(inc_chain == botVertex) {
1046     reflexChain rChain(20, 0);
1047     rChain.processNewVertex(topVertex->getVertex(top_index), pStream);
1048     for(i=dec_index; i< dec_chain->get_npoints(); i++){
1049       rChain.processNewVertex(dec_chain->getVertex(i), pStream);
1050     }
1051     for(temp = dec_chain->getPrev(); temp != botVertex; temp = temp->getPrev())
1052       {
1053         for(i=0; i<temp->get_npoints(); i++){
1054           rChain.processNewVertex(temp->getVertex(i), pStream);
1055         }
1056       }
1057   }
1058   else if(dec_chain==botVertex) {
1059     reflexChain rChain(20, 1);
1060     rChain.processNewVertex(topVertex->getVertex(top_index), pStream);
1061     for(i=inc_index; i< inc_chain->get_npoints(); i++){
1062       rChain.processNewVertex(inc_chain->getVertex(i), pStream);
1063     }
1064     for(temp = inc_chain->getPrev(); temp != botVertex; temp = temp->getNext())
1065       {
1066         for(i=0; i<temp->get_npoints(); i++){
1067           rChain.processNewVertex(temp->getVertex(i), pStream);
1068         }
1069       }
1070   }
1071   else /*neither reached the bottom*/{
1072     if(compV2InY(inc_chain->getVertex(inc_index), dec_chain->getVertex(dec_index)) <=0) {
1073       reflexChain rChain(20, 0);
1074       rChain.processNewVertex(topVertex -> getVertex(top_index), pStream);
1075       temp = dec_chain;
1076       tempIndex = dec_index;
1077       while( compV2InY(inc_chain->getVertex(inc_index), temp->getVertex(tempIndex))<=0) {
1078         oldtemp = temp;
1079         oldtempIndex = tempIndex;
1080         rChain.processNewVertex(temp->getVertex(tempIndex), pStream);
1081
1082         if(tempIndex == temp->get_npoints()-1){
1083           tempIndex = 0;
1084           temp = temp->getPrev();
1085         }
1086         else{
1087           tempIndex++;
1088         }
1089       }
1090       rChain.outputFan(inc_chain->getVertex(inc_index), pStream);
1091       monoTriangulationRec(inc_chain, inc_index, temp, tempIndex, oldtemp, oldtempIndex, botVertex, pStream);
1092     }
1093     else /* >0*/ {
1094       reflexChain rChain(20, 1);
1095       rChain.processNewVertex(topVertex -> getVertex(top_index), pStream);
1096       temp = inc_chain;
1097       tempIndex = inc_index;
1098       while( compV2InY(temp->getVertex(tempIndex), dec_chain->getVertex(dec_index))>0){
1099         oldtemp = temp;
1100         oldtempIndex = tempIndex;
1101         rChain.processNewVertex(temp->getVertex(tempIndex), pStream);
1102
1103         if(tempIndex == temp->get_npoints()-1){
1104           tempIndex = 0;
1105           temp = temp->getNext();
1106         }
1107         else{
1108           tempIndex++;
1109         }
1110       }
1111       rChain.outputFan(dec_chain->getVertex(dec_index), pStream);
1112       monoTriangulationRec(temp, tempIndex, dec_chain, dec_index, oldtemp, oldtempIndex, botVertex, pStream);
1113     }
1114   } /*end case neither reached the bottom*/
1115 }
1116
1117 /***************************vertexArray begin here**********************************/
1118 vertexArray::vertexArray(Real2* vertices, Int nVertices)
1119 {
1120   Int i;
1121   size = index = nVertices;
1122   array = (Real**) malloc(sizeof(Real*) * nVertices);
1123   assert(array);
1124   for(i=0; i<nVertices; i++)
1125     {
1126       array[i] = vertices[i];
1127       array[i] = vertices[i];
1128     }
1129 }
1130
1131 vertexArray::vertexArray(Int s)
1132 {
1133   size = s;
1134   array = (Real**) malloc(sizeof(Real*) * s);
1135   assert(array);
1136   index = 0;
1137 }
1138
1139 vertexArray::~vertexArray()
1140 {
1141   free(array);
1142 }
1143
1144 void vertexArray::appendVertex(Real* ptr)
1145 {
1146   Int i;
1147   if(index >= size){
1148     Real** temp = (Real**) malloc(sizeof(Real*) * (2*size +1));
1149     assert(temp);
1150     for(i=0; i<index; i++)
1151       temp[i] = array[i];
1152     free(array);
1153     array = temp;
1154     size = 2*size+1;
1155   }
1156   array[index++] = ptr;
1157 }
1158
1159 void vertexArray::print()
1160 {
1161   printf("vertex Array:index=%i, size=%i\n", index, size);
1162   for(Int i=0; i<index; i++)
1163     {
1164       printf("(%f,%f) ", array[i][0], array[i][1]);
1165     }
1166   printf("\n");
1167 }
1168
1169 /*find the first i such that array[i][1] >= v
1170  * and array[i+1][1] <v
1171  * if index == 0 (the array is empty, return -1.
1172  * if v is above all, return -1.
1173  * if v is below all, return index-1.
1174  */
1175 Int vertexArray::findIndexAbove(Real v)
1176 {
1177   Int i;
1178   if(index == 0)
1179     return -1;
1180   else if(array[0][1] < v)
1181     return -1;
1182   else
1183     {
1184       for(i=1; i<index; i++)
1185         {
1186           if(array[i][1] < v)
1187             break;
1188         }
1189       return i-1;
1190     }
1191 }
1192
1193 /*find the first i<=endIndex such that array[i][1] <= v
1194  * and array[i-1][1] > v
1195  *if sartIndex>endIndex, then return endIndex+1.
1196  *otherwise, startIndex<=endIndex, it is assumed that
1197  * 0<=startIndex<=endIndex<index.
1198  * if v is below all, return endIndex+1
1199  * if v is above all, return startIndex.
1200  */
1201 Int vertexArray::findIndexBelowGen(Real v, Int startIndex, Int endIndex)
1202 {
1203   Int i;
1204   if(startIndex > endIndex)
1205     return endIndex+1;
1206   else if(array[endIndex][1] >  v)
1207     return endIndex+1;
1208   else //now array[endIndex][1] <= v
1209     {
1210       for(i=endIndex-1; i>=startIndex; i--)
1211         {
1212           if(array[i][1] > v)
1213             break;
1214         }
1215       return i+1;
1216     }
1217 }
1218
1219 /*find the first i<=endIndex such that array[i-1][1] >= v
1220  * and array[i][1] < v
1221  *if sartIndex>endIndex, then return endIndex+1.
1222  *otherwise, startIndex<=endIndex, it is assumed that
1223  * 0<=startIndex<=endIndex<index.
1224  * if v is below or equal to all, return endIndex+1
1225  * if v is strictly above all, return startIndex.
1226  */
1227 Int vertexArray::findIndexStrictBelowGen(Real v, Int startIndex, Int endIndex)
1228 {
1229   Int i;
1230   if(startIndex > endIndex)
1231     return endIndex+1;
1232   else if(array[endIndex][1] >=  v)
1233     return endIndex+1;
1234   else //now array[endIndex][1] < v
1235     {
1236       for(i=endIndex-1; i>=startIndex; i--)
1237         {
1238           if(array[i][1] >= v)
1239             break;
1240         }
1241       return i+1;
1242     }
1243 }
1244
1245 /*find the first i>startIndex such that array[i-1][1] > v
1246  * and array[i][1] >=v
1247  *if sartIndex>endIndex, then return startIndex-1.
1248  *otherwise, startIndex<=endIndex, it is assumed that
1249  * 0<=startIndex<=endIndex<index.
1250  * if v is strictly above all, return startIndex-1
1251  * if v is strictly below all, return endIndex.
1252  */
1253 Int vertexArray::findIndexFirstAboveEqualGen(Real v, Int startIndex, Int endIndex)
1254 {
1255
1256   Int i;
1257   if(startIndex > endIndex)
1258     return startIndex-1;
1259   else if(array[startIndex][1] < v)
1260     return startIndex-1;
1261   else //now array[startIndex][1] >= v
1262     {
1263
1264       for(i=startIndex; i<=endIndex; i++)
1265         {
1266           if(array[i][1] <= v)
1267             break;
1268         }
1269       if(i>endIndex) // v is strictly below all
1270         return endIndex;
1271       else if(array[i][1] == v)
1272         return i;
1273       else
1274         return i-1;
1275     }
1276
1277 }
1278
1279
1280 /*find the first i>=startIndex such that array[i][1] >= v
1281  * and array[i+1][1] <v
1282  *if sartIndex>endIndex, then return startIndex-1.
1283  *otherwise, startIndex<=endIndex, it is assumed that
1284  * 0<=startIndex<=endIndex<index.
1285  * if v is above all, return startIndex-1
1286  * if v is below all, return endIndex.
1287  */
1288 Int vertexArray::findIndexAboveGen(Real v, Int startIndex, Int endIndex)
1289 {
1290   Int i;
1291   if(startIndex > endIndex)
1292     return startIndex-1;
1293   else if(array[startIndex][1] < v)
1294     return startIndex-1;
1295   else //now array[startIndex][1] >= v
1296     {
1297       for(i=startIndex+1; i<=endIndex; i++)
1298         {
1299           if(array[i][1] < v)
1300             break;
1301         }
1302       return i-1;
1303     }
1304 }
1305
1306 Int vertexArray::findDecreaseChainFromEnd(Int begin, Int end)
1307 {
1308   Int i = end;
1309   Real prevU = array[i][0];
1310   Real thisU;
1311   for(i=end-1; i>=begin; i--){
1312     thisU = array[i][0];
1313     if(thisU < prevU)
1314       prevU = thisU;
1315     else
1316       break;
1317   }
1318   return i;
1319 }
1320
1321 //if(V(start) == v, return start, other wise return the
1322 //last i so that V(i)==v
1323 Int vertexArray::skipEqualityFromStart(Real v, Int start, Int end)
1324 {
1325   Int i;
1326   if(array[start][1] != v)
1327     return start;
1328   //now array[start][1] == v
1329   for(i=start+1; i<= end; i++)
1330     if(array[i][1] != v)
1331       break;
1332   return i-1;
1333 }
1334
1335
1336 /***************************vertexArray end****************************************/
1337
1338
1339
1340 /***************************relfex chain stuff begin here*****************************/
1341
1342 reflexChain::reflexChain(Int size, Int is_increasing)
1343 {
1344   queue = (Real2*) malloc(sizeof(Real2) * size);
1345   assert(queue);
1346   index_queue = 0;
1347   size_queue = size;
1348   isIncreasing = is_increasing;
1349 }
1350
1351 reflexChain::~reflexChain()
1352 {
1353   free(queue);
1354 }
1355
1356 /*put (u,v) at the end of the queue
1357  *pay attention to space
1358  */
1359 void reflexChain::insert(Real u, Real v)
1360 {
1361   Int i;
1362   if(index_queue >= size_queue) {
1363     Real2 *temp = (Real2*) malloc(sizeof(Real2) * (2*size_queue+1));
1364     assert(temp);
1365
1366     /*copy*/
1367     for(i=0; i<index_queue; i++){
1368       temp[i][0] = queue[i][0];
1369       temp[i][1] = queue[i][1];
1370     }
1371
1372     free(queue);
1373     queue = temp;
1374     size_queue = 2*size_queue + 1;
1375   }
1376
1377   queue[index_queue][0] = u;
1378   queue[index_queue][1] = v;
1379   index_queue ++;
1380 }
1381
1382 void reflexChain::insert(Real v[2])
1383 {
1384   insert(v[0], v[1]);
1385 }
1386
1387 /*
1388 static Real area(Real A[2], Real B[2], Real C[2])
1389 {
1390   Real Bx, By, Cx, Cy;
1391   Bx = B[0] - A[0];
1392   By = B[1] - A[1];
1393   Cx = C[0] - A[0];
1394   Cy = C[1] - A[1];
1395   return Bx*Cy - Cx*By;
1396 }
1397 */
1398
1399 /*the chain is reflex, and the vertex v is
1400  *on the other side of the chain, so that
1401  *we can outout the fan with v as the
1402  *the center
1403  */
1404 void reflexChain::outputFan(Real v[2], primStream* pStream)
1405 {
1406   Int i;
1407   pStream->begin();
1408   pStream->insert(v);
1409   if(isIncreasing) {
1410     for(i=0; i<index_queue; i++)
1411       pStream->insert(queue[i]);
1412   }
1413   else {
1414     for(i=index_queue-1; i>=0; i--)
1415       pStream->insert(queue[i]);
1416   }
1417   pStream->end(PRIMITIVE_STREAM_FAN);
1418 }
1419
1420 void reflexChain::processNewVertex(Real v[2], primStream* pStream)
1421 {
1422   Int i,j,k;
1423   Int isReflex;
1424   /*if there are at most one vertex in the queue, then simply insert
1425    */
1426   if(index_queue <=1){
1427     insert(v);
1428     return;
1429   }
1430
1431   /*there are at least two vertices in the queue*/
1432   j=index_queue-1;
1433
1434   for(i=j; i>=1; i--) {
1435     if(isIncreasing) {
1436       isReflex = (area(queue[i-1], queue[i], v) <= 0.0);
1437     }
1438     else /*decreasing*/{
1439       isReflex = (area(v, queue[i], queue[i-1]) <= 0.0);
1440     }
1441     if(isReflex) {
1442       break;
1443     }
1444   }
1445
1446   /*
1447    *if i<j then vertices: i+1--j are convex
1448    * output triangle fan:
1449    *  v, and queue[i], i+1, ..., j
1450    */
1451   if(i<j)
1452     {
1453       pStream->begin();
1454       pStream->insert(v);
1455       if(isIncreasing) {
1456         for(k=i; k<=j; k++)
1457           pStream->insert(queue[k]);
1458       }
1459       else {
1460         for(k=j; k>=i; k--)
1461           pStream->insert(queue[k]);
1462       }
1463
1464       pStream->end(PRIMITIVE_STREAM_FAN);
1465     }
1466
1467   /*delete vertices i+1--j from the queue*/
1468   index_queue = i+1;
1469   /*finally insert v at the end of the queue*/
1470   insert(v);
1471
1472 }
1473
1474 void reflexChain::print()
1475 {
1476   Int i;
1477   printf("reflex chain: isIncreasing=%i\n", isIncreasing);
1478   for(i=0; i<index_queue; i++) {
1479     printf("(%f,%f) ", queue[i][0], queue[i][1]);
1480   }
1481   printf("\n");
1482 }