src/glu/sgi/libnurbs/nurbtess/sampleMonoPoly.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 <math.h>
  42
  43 #ifndef max
  44 #define max(a,b) ((a>b)? a:b)
  45 #endif
  46 #ifndef min
  47 #define min(a,b) ((a>b)? b:a)
  48 #endif
  49
  50 #include <GL/gl.h>
  51
  52 #include "glimports.h"
  53 #include "zlassert.h"
  54 #include "sampleMonoPoly.h"
  55 #include "sampleComp.h"
  56 #include "polyDBG.h"
  57 #include "partitionX.h"
  58
  59
  60 #define ZERO 0.00001
  61
  62 //#define  MYDEBUG
  63
  64 //#define SHORTEN_GRID_LINE
  65 //see work/newtess/internal/test/problems
  66
  67
  68 /*split a polygon so that each vertex correcpond to one edge
  69  *the head of the first edge of the returned plygon must be the head of the first
  70  *edge of the origianl polygon. This is crucial for the code in sampleMonoPoly function
  71  */
  72  directedLine*  polygonConvert(directedLine* polygon)
  73 {
  74   int i;
  75   directedLine* ret;
  76   sampledLine* sline;
  77   sline = new sampledLine(2);
  78   sline->setPoint(0, polygon->getVertex(0));
  79   sline->setPoint(1, polygon->getVertex(1));
  80   ret=new directedLine(INCREASING, sline);
  81   for(i=1; i<= polygon->get_npoints()-2; i++)
  82     {
  83       sline = new sampledLine(2);
  84       sline->setPoint(0, polygon->getVertex(i));
  85       sline->setPoint(1, polygon->getVertex(i+1));
  86       ret->insert(new directedLine(INCREASING, sline));
  87     }
  88
  89   for(directedLine *temp = polygon->getNext(); temp != polygon; temp = temp->getNext())
  90     {
  91       for(i=0; i<= temp->get_npoints()-2; i++)
  92         {
  93           sline = new sampledLine(2);
  94           sline->setPoint(0, temp->getVertex(i));
  95           sline->setPoint(1, temp->getVertex(i+1));
  96           ret->insert(new directedLine(INCREASING, sline));
  97         }
  98     }
  99   return ret;
 100 }
 101
 102 void triangulateConvexPolyVertical(directedLine* topV, directedLine* botV, primStream *pStream)
 103 {
 104   Int i,j;
 105   Int n_leftVerts;
 106   Int n_rightVerts;
 107   Real** leftVerts;
 108   Real** rightVerts;
 109   directedLine* tempV;
 110   n_leftVerts = 0;
 111   for(tempV = topV; tempV != botV; tempV = tempV->getNext())
 112     {
 113       n_leftVerts += tempV->get_npoints();
 114     }
 115   n_rightVerts=0;
 116   for(tempV = botV; tempV != topV; tempV = tempV->getNext())
 117     {
 118       n_rightVerts += tempV->get_npoints();
 119     }
 120
 121   Real2* temp_leftVerts = (Real2 *) malloc(sizeof(Real2) * n_leftVerts);
 122   assert(temp_leftVerts);
 123   Real2* temp_rightVerts = (Real2 *) malloc(sizeof(Real2) * n_rightVerts);
 124   assert(temp_rightVerts);
 125
 126   leftVerts = (Real**) malloc(sizeof(Real2*) * n_leftVerts);
 127   assert(leftVerts);
 128   rightVerts = (Real**) malloc(sizeof(Real2*) * n_rightVerts);
 129   assert(rightVerts);
 130   for(i=0; i<n_leftVerts; i++)
 131     leftVerts[i] = temp_leftVerts[i];
 132   for(i=0; i<n_rightVerts; i++)
 133     rightVerts[i] = temp_rightVerts[i];
 134
 135   i=0;
 136   for(tempV = topV; tempV != botV; tempV = tempV->getNext())
 137     {
 138       for(j=1; j<tempV->get_npoints(); j++)
 139         {
 140           leftVerts[i][0] = tempV->getVertex(j)[0];
 141           leftVerts[i][1] = tempV->getVertex(j)[1];
 142           i++;
 143         }
 144     }
 145   n_leftVerts = i;
 146   i=0;
 147   for(tempV = topV->getPrev(); tempV != botV->getPrev(); tempV = tempV->getPrev())
 148     {
 149       for(j=tempV->get_npoints()-1; j>=1; j--)
 150         {
 151           rightVerts[i][0] = tempV->getVertex(j)[0];
 152           rightVerts[i][1] = tempV->getVertex(j)[1];
 153           i++;
 154         }
 155     }
 156   n_rightVerts = i;
 157   triangulateXYMonoTB(n_leftVerts, leftVerts, n_rightVerts, rightVerts, pStream);
 158   free(leftVerts);
 159   free(rightVerts);
 160   free(temp_leftVerts);
 161   free(temp_rightVerts);
 162 }
 163
 164 void triangulateConvexPolyHoriz(directedLine* leftV, directedLine* rightV, primStream *pStream)
 165 {
 166   Int i,j;
 167   Int n_lowerVerts;
 168   Int n_upperVerts;
 169   Real2 *lowerVerts;
 170   Real2 *upperVerts;
 171   directedLine* tempV;
 172   n_lowerVerts=0;
 173   for(tempV = leftV; tempV != rightV; tempV = tempV->getNext())
 174     {
 175       n_lowerVerts += tempV->get_npoints();
 176     }
 177   n_upperVerts=0;
 178   for(tempV = rightV; tempV != leftV; tempV = tempV->getNext())
 179     {
 180       n_upperVerts += tempV->get_npoints();
 181     }
 182   lowerVerts = (Real2 *) malloc(sizeof(Real2) * n_lowerVerts);
 183   assert(n_lowerVerts);
 184   upperVerts = (Real2 *) malloc(sizeof(Real2) * n_upperVerts);
 185   assert(n_upperVerts);
 186   i=0;
 187   for(tempV = leftV; tempV != rightV; tempV = tempV->getNext())
 188     {
 189       for(j=0; j<tempV->get_npoints(); j++)
 190         {
 191           lowerVerts[i][0] = tempV->getVertex(j)[0];
 192           lowerVerts[i][1] = tempV->getVertex(j)[1];
 193           i++;
 194         }
 195     }
 196   i=0;
 197   for(tempV = leftV->getPrev(); tempV != rightV->getPrev(); tempV = tempV->getPrev())
 198     {
 199       for(j=tempV->get_npoints()-1; j>=0; j--)
 200         {
 201           upperVerts[i][0] = tempV->getVertex(j)[0];
 202           upperVerts[i][1] = tempV->getVertex(j)[1];
 203           i++;
 204         }
 205     }
 206   triangulateXYMono(n_upperVerts, upperVerts, n_lowerVerts, lowerVerts, pStream);
 207   free(lowerVerts);
 208   free(upperVerts);
 209 }
 210 void triangulateConvexPoly(directedLine* polygon, Int ulinear, Int vlinear, primStream* pStream)
 211 {
 212   /*find left, right, top , bot
 213     */
 214   directedLine* tempV;
 215   directedLine* topV;
 216   directedLine* botV;
 217   directedLine* leftV;
 218   directedLine* rightV;
 219   topV = botV = polygon;
 220
 221   for(tempV = polygon->getNext(); tempV != polygon; tempV = tempV->getNext())
 222     {
 223       if(compV2InY(topV->head(), tempV->head())<0) {
 224
 225         topV = tempV;
 226       }
 227       if(compV2InY(botV->head(), tempV->head())>0) {
 228
 229         botV = tempV;
 230       }
 231     }
 232   //find leftV
 233   for(tempV = topV; tempV != botV; tempV = tempV->getNext())
 234     {
 235       if(tempV->tail()[0] >= tempV->head()[0])
 236         break;
 237     }
 238   leftV = tempV;
 239   //find rightV
 240   for(tempV = botV; tempV != topV; tempV = tempV->getNext())
 241     {
 242       if(tempV->tail()[0] <= tempV->head()[0])
 243         break;
 244     }
 245   rightV = tempV;
 246   if(vlinear)
 247     {
 248       triangulateConvexPolyHoriz( leftV, rightV, pStream);
 249     }
 250   else if(ulinear)
 251     {
 252       triangulateConvexPolyVertical(topV, botV, pStream);
 253     }
 254   else
 255     {
 256       if(DBG_is_U_direction(polygon))
 257         {
 258           triangulateConvexPolyHoriz( leftV, rightV, pStream);
 259         }
 260       else
 261         triangulateConvexPolyVertical(topV, botV, pStream);
 262     }
 263 }
 264
 265 /*for debug purpose*/
 266 void drawCorners(
 267                  Real* topV, Real* botV,
 268                  vertexArray* leftChain,
 269                  vertexArray* rightChain,
 270                  gridBoundaryChain* leftGridChain,
 271                  gridBoundaryChain* rightGridChain,
 272                  Int gridIndex1,
 273                  Int gridIndex2,
 274                  Int leftCornerWhere,
 275                  Int leftCornerIndex,
 276                  Int rightCornerWhere,
 277                  Int rightCornerIndex,
 278                  Int bot_leftCornerWhere,
 279                  Int bot_leftCornerIndex,
 280                  Int bot_rightCornerWhere,
 281                  Int bot_rightCornerIndex)
 282 {
 283   Real* leftCornerV;
 284   Real* rightCornerV;
 285   Real* bot_leftCornerV;
 286   Real* bot_rightCornerV;
 287
 288   if(leftCornerWhere == 1)
 289     leftCornerV = topV;
 290   else if(leftCornerWhere == 0)
 291     leftCornerV = leftChain->getVertex(leftCornerIndex);
 292   else
 293     leftCornerV = rightChain->getVertex(leftCornerIndex);
 294
 295   if(rightCornerWhere == 1)
 296     rightCornerV = topV;
 297   else if(rightCornerWhere == 0)
 298     rightCornerV = leftChain->getVertex(rightCornerIndex);
 299   else
 300     rightCornerV = rightChain->getVertex(rightCornerIndex);
 301
 302   if(bot_leftCornerWhere == 1)
 303     bot_leftCornerV = botV;
 304   else if(bot_leftCornerWhere == 0)
 305     bot_leftCornerV = leftChain->getVertex(bot_leftCornerIndex);
 306   else
 307     bot_leftCornerV = rightChain->getVertex(bot_leftCornerIndex);
 308
 309   if(bot_rightCornerWhere == 1)
 310     bot_rightCornerV = botV;
 311   else if(bot_rightCornerWhere == 0)
 312     bot_rightCornerV = leftChain->getVertex(bot_rightCornerIndex);
 313   else
 314     bot_rightCornerV = rightChain->getVertex(bot_rightCornerIndex);
 315
 316   Real topGridV = leftGridChain->get_v_value(gridIndex1);
 317   Real topGridU1 = leftGridChain->get_u_value(gridIndex1);
 318   Real topGridU2 = rightGridChain->get_u_value(gridIndex1);
 319   Real botGridV = leftGridChain->get_v_value(gridIndex2);
 320   Real botGridU1 = leftGridChain->get_u_value(gridIndex2);
 321   Real botGridU2 = rightGridChain->get_u_value(gridIndex2);
 322
 323   glBegin(GL_LINE_STRIP);
 324   glVertex2fv(leftCornerV);
 325   glVertex2f(topGridU1, topGridV);
 326   glEnd();
 327
 328   glBegin(GL_LINE_STRIP);
 329   glVertex2fv(rightCornerV);
 330   glVertex2f(topGridU2, topGridV);
 331   glEnd();
 332
 333   glBegin(GL_LINE_STRIP);
 334   glVertex2fv(bot_leftCornerV);
 335   glVertex2f(botGridU1, botGridV);
 336   glEnd();
 337
 338   glBegin(GL_LINE_STRIP);
 339   glVertex2fv(bot_rightCornerV);
 340   glVertex2f(botGridU2, botGridV);
 341   glEnd();
 342
 343
 344 }
 345
 346 void toVertexArrays(directedLine* topV, directedLine* botV, vertexArray& leftChain, vertexArray& rightChain)
 347 {
 348   Int i;
 349   directedLine* tempV;
 350   for(i=1; i<=topV->get_npoints()-2; i++) { /*the first vertex is the top vertex which doesn't belong to inc_chain*/
 351     leftChain.appendVertex(topV->getVertex(i));
 352   }
 353   for(tempV = topV->getNext(); tempV != botV; tempV = tempV->getNext())
 354     {
 355       for(i=0; i<=tempV->get_npoints()-2; i++){
 356         leftChain.appendVertex(tempV->getVertex(i));
 357       }
 358     }
 359
 360   for(tempV = topV->getPrev(); tempV != botV; tempV = tempV->getPrev())
 361     {
 362       for(i=tempV->get_npoints()-2; i>=0; i--){
 363         rightChain.appendVertex(tempV->getVertex(i));
 364       }
 365     }
 366   for(i=botV->get_npoints()-2; i>=1; i--){
 367     rightChain.appendVertex(tempV->getVertex(i));
 368   }
 369 }
 370
 371
 372 void findTopAndBot(directedLine* polygon, directedLine*& topV, directedLine*& botV)
 373 {
 374   assert(polygon);
 375   directedLine* tempV;
 376   topV = botV = polygon;
 377    for(tempV = polygon->getNext(); tempV != polygon; tempV = tempV->getNext())
 378     {
 379       if(compV2InY(topV->head(), tempV->head())<0) {
 380         topV = tempV;
 381       }
 382       if(compV2InY(botV->head(), tempV->head())>0) {
 383         botV = tempV;
 384       }
 385     }
 386 }
 387
 388 void findGridChains(directedLine* topV, directedLine* botV,
 389                     gridWrap* grid,
 390                     gridBoundaryChain*& leftGridChain,
 391                     gridBoundaryChain*& rightGridChain)
 392 {
 393   /*find the first(top) and the last (bottom) grid line which intersect the
 394    *this polygon
 395    */
 396   Int firstGridIndex; /*the index in the grid*/
 397   Int lastGridIndex;
 398
 399   firstGridIndex = (Int) ((topV->head()[1] - grid->get_v_min()) / (grid->get_v_max() - grid->get_v_min()) * (grid->get_n_vlines()-1));
 400
 401   if(botV->head()[1] < grid->get_v_min())
 402     lastGridIndex = 0;
 403   else
 404     lastGridIndex  = (Int) ((botV->head()[1] - grid->get_v_min()) / (grid->get_v_max() - grid->get_v_min()) * (grid->get_n_vlines()-1)) + 1;
 405
 406   /*find the interval inside  the polygon for each gridline*/
 407   Int *leftGridIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
 408   assert(leftGridIndices);
 409   Int *rightGridIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
 410   assert(rightGridIndices);
 411   Int *leftGridInnerIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
 412   assert(leftGridInnerIndices);
 413   Int *rightGridInnerIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
 414   assert(rightGridInnerIndices);
 415
 416   findLeftGridIndices(topV, firstGridIndex, lastGridIndex, grid,  leftGridIndices, leftGridInnerIndices);
 417
 418   findRightGridIndices(topV, firstGridIndex, lastGridIndex, grid,  rightGridIndices, rightGridInnerIndices);
 419
 420   leftGridChain =  new gridBoundaryChain(grid, firstGridIndex, firstGridIndex-lastGridIndex+1, leftGridIndices, leftGridInnerIndices);
 421
 422   rightGridChain = new gridBoundaryChain(grid, firstGridIndex, firstGridIndex-lastGridIndex+1, rightGridIndices, rightGridInnerIndices);
 423
 424   free(leftGridIndices);
 425   free(rightGridIndices);
 426   free(leftGridInnerIndices);
 427   free(rightGridInnerIndices);
 428 }
 429
 430 void findDownCorners(Real *botVertex,
 431                    vertexArray *leftChain, Int leftChainStartIndex, Int leftChainEndIndex,
 432                    vertexArray *rightChain, Int rightChainStartIndex, Int rightChainEndIndex,
 433                    Real v,
 434                    Real uleft,
 435                    Real uright,
 436                    Int& ret_leftCornerWhere, /*0: left chain, 1: topvertex, 2: rightchain*/
 437                    Int& ret_leftCornerIndex, /*useful when ret_leftCornerWhere == 0 or 2*/
 438                    Int& ret_rightCornerWhere, /*0: left chain, 1: topvertex, 2: rightchain*/
 439                    Int& ret_rightCornerIndex /*useful when ret_leftCornerWhere == 0 or 2*/
 440                    )
 441 {
 442 #ifdef MYDEBUG
 443 printf("*************enter find donw corner\n");
 444 printf("finddownCorner: v=%f, uleft=%f, uright=%f\n", v, uleft, uright);
 445 printf("(%i,%i,%i,%i)\n", leftChainStartIndex, leftChainEndIndex,rightChainStartIndex, rightChainEndIndex);
 446 printf("left chain is\n");
 447 leftChain->print();
 448 printf("right chain is\n");
 449 rightChain->print();
 450 #endif
 451
 452   assert(v > botVertex[1]);
 453   Real leftGridPoint[2];
 454   leftGridPoint[0] = uleft;
 455   leftGridPoint[1] = v;
 456   Real rightGridPoint[2];
 457   rightGridPoint[0] = uright;
 458   rightGridPoint[1] = v;
 459
 460   Int i;
 461   Int index1, index2;
 462
 463   index1 = leftChain->findIndexBelowGen(v, leftChainStartIndex, leftChainEndIndex);
 464   index2 = rightChain->findIndexBelowGen(v, rightChainStartIndex, rightChainEndIndex);
 465
 466   if(index2 <= rightChainEndIndex) //index2 was found above
 467     index2 = rightChain->skipEqualityFromStart(v, index2, rightChainEndIndex);
 468
 469   if(index1>leftChainEndIndex && index2 > rightChainEndIndex) /*no point below v on left chain or right chain*/
 470     {
 471
 472       /*the botVertex is the only vertex below v*/
 473       ret_leftCornerWhere = 1;
 474       ret_rightCornerWhere = 1;
 475     }
 476   else if(index1>leftChainEndIndex ) /*index2 <= rightChainEndIndex*/
 477     {
 478
 479       ret_rightCornerWhere = 2; /*on right chain*/
 480       ret_rightCornerIndex = index2;
 481
 482
 483       Real tempMin = rightChain->getVertex(index2)[0];
 484       Int tempI = index2;
 485       for(i=index2+1; i<= rightChainEndIndex; i++)
 486         if(rightChain->getVertex(i)[0] < tempMin)
 487           {
 488             tempI = i;
 489             tempMin = rightChain->getVertex(i)[0];
 490           }
 491
 492
 493       //we consider whether we can use botVertex as left corner. First check
 494       //if (leftGirdPoint, botVertex) interesects right chian or not.
 495      if(DBG_intersectChain(rightChain, rightChainStartIndex,rightChainEndIndex,
 496                                     leftGridPoint, botVertex))
 497        {
 498          ret_leftCornerWhere = 2;//right
 499          ret_leftCornerIndex = index2; //should use tempI???
 500        }
 501      else if(botVertex[0] < tempMin)
 502        ret_leftCornerWhere = 1; //bot
 503      else
 504        {
 505          ret_leftCornerWhere = 2; //right
 506          ret_leftCornerIndex = tempI;
 507        }
 508     }
 509   else if(index2> rightChainEndIndex) /*index1<=leftChainEndIndex*/
 510     {
 511       ret_leftCornerWhere = 0; /*left chain*/
 512       ret_leftCornerIndex = index1;
 513
 514       /*find the vertex on the left chain with the maximum u,
 515        *either this vertex or the botvertex can be used as the right corner
 516        */
 517
 518       Int tempI;
 519       //skip those points which are equal to v. (avoid degeneratcy)
 520       for(tempI = index1; tempI <= leftChainEndIndex; tempI++)
 521         if(leftChain->getVertex(tempI)[1] < v)
 522           break;
 523       if(tempI > leftChainEndIndex)
 524         ret_rightCornerWhere = 1;
 525       else
 526         {
 527           Real tempMax = leftChain->getVertex(tempI)[0];
 528           for(i=tempI; i<= leftChainEndIndex; i++)
 529             if(leftChain->getVertex(i)[0] > tempMax)
 530               {
 531                 tempI = i;
 532                 tempMax = leftChain->getVertex(i)[0];
 533               }
 534
 535
 536
 537           //we consider whether we can use botVertex as a corner. So first we check
 538           //whether (rightGridPoint, botVertex) interescts the left chain or not.
 539           if(DBG_intersectChain(leftChain, leftChainStartIndex,leftChainEndIndex,
 540                                     rightGridPoint, botVertex))
 541             {
 542               ret_rightCornerWhere = 0;
 543               ret_rightCornerIndex = index1; //should use tempI???
 544             }
 545           else if(botVertex[0] > tempMax)
 546             {
 547
 548               ret_rightCornerWhere = 1;
 549             }
 550           else
 551             {
 552               ret_rightCornerWhere = 0;
 553               ret_rightCornerIndex = tempI;
 554             }
 555         }
 556
 557     }
 558   else /*index1<=leftChainEndIndex and index2 <=rightChainEndIndex*/
 559     {
 560       if(leftChain->getVertex(index1)[1] >= rightChain->getVertex(index2)[1]) /*left point above right point*/
 561         {
 562           ret_leftCornerWhere = 0; /*on left chain*/
 563           ret_leftCornerIndex = index1;
 564
 565           Real tempMax;
 566           Int tempI;
 567
 568           tempI = index1;
 569           tempMax = leftChain->getVertex(index1)[0];
 570
 571           /*find the maximum u for all the points on the left above the right point index2*/
 572           for(i=index1+1; i<= leftChainEndIndex; i++)
 573             {
 574               if(leftChain->getVertex(i)[1] < rightChain->getVertex(index2)[1])
 575                 break;
 576
 577               if(leftChain->getVertex(i)[0]>tempMax)
 578                 {
 579                   tempI = i;
 580                   tempMax = leftChain->getVertex(i)[0];
 581                 }
 582             }
 583           //we consider if we can use rightChain(index2) as right corner
 584           //we check if (rightChain(index2), rightGidPoint) intersecs left chain or not.
 585           if(DBG_intersectChain(leftChain, leftChainStartIndex,leftChainEndIndex, rightGridPoint, rightChain->getVertex(index2)))
 586             {
 587               ret_rightCornerWhere = 0;
 588               ret_rightCornerIndex = index1; //should use tempI???
 589             }
 590           else if(tempMax >= rightChain->getVertex(index2)[0] ||
 591              tempMax >= uright
 592              )
 593             {
 594
 595               ret_rightCornerWhere = 0; /*on left Chain*/
 596               ret_rightCornerIndex = tempI;
 597             }
 598           else
 599             {
 600               ret_rightCornerWhere = 2; /*on right chain*/
 601               ret_rightCornerIndex = index2;
 602             }
 603         }
 604       else /*left below right*/
 605         {
 606           ret_rightCornerWhere = 2; /*on the right*/
 607           ret_rightCornerIndex = index2;
 608
 609           Real tempMin;
 610           Int tempI;
 611
 612           tempI = index2;
 613           tempMin = rightChain->getVertex(index2)[0];
 614
 615           /*find the minimum u for all the points on the right above the left poitn index1*/
 616           for(i=index2+1; i<= rightChainEndIndex; i++)
 617             {
 618               if( rightChain->getVertex(i)[1] < leftChain->getVertex(index1)[1])
 619                 break;
 620               if(rightChain->getVertex(i)[0] < tempMin)
 621                 {
 622                   tempI = i;
 623                   tempMin = rightChain->getVertex(i)[0];
 624                 }
 625             }
 626
 627           //we consider if we can use leftchain(index1) as left corner.
 628           //we check if (leftChain(index1) intersects right chian or not
 629           if(DBG_intersectChain(rightChain, rightChainStartIndex, rightChainEndIndex, leftGridPoint, leftChain->getVertex(index1)))
 630             {
 631               ret_leftCornerWhere = 2;
 632               ret_leftCornerIndex = index2; //should use tempI???
 633               }
 634           else if(tempMin <= leftChain->getVertex(index1)[0] ||
 635              tempMin <= uleft)
 636             {
 637               ret_leftCornerWhere = 2; /* on right chain*/
 638               ret_leftCornerIndex = tempI;
 639             }
 640           else
 641             {
 642               ret_leftCornerWhere = 0; /*on left chain*/
 643               ret_leftCornerIndex = index1;
 644             }
 645         }
 646     }
 647
 648 }
 649
 650
 651 void findUpCorners(Real *topVertex,
 652                    vertexArray *leftChain, Int leftChainStartIndex, Int leftChainEndIndex,
 653                    vertexArray *rightChain, Int rightChainStartIndex, Int rightChainEndIndex,
 654                    Real v,
 655                    Real uleft,
 656                    Real uright,
 657                    Int& ret_leftCornerWhere, /*0: left chain, 1: topvertex, 2: rightchain*/
 658                    Int& ret_leftCornerIndex, /*useful when ret_leftCornerWhere == 0 or 2*/
 659                    Int& ret_rightCornerWhere, /*0: left chain, 1: topvertex, 2: rightchain*/
 660                    Int& ret_rightCornerIndex /*useful when ret_leftCornerWhere == 0 or 2*/
 661                    )
 662 {
 663 #ifdef MYDEBUG
 664 printf("***********enter findUpCorners\n");
 665 #endif
 666
 667   assert(v < topVertex[1]);
 668   Real leftGridPoint[2];
 669   leftGridPoint[0] = uleft;
 670   leftGridPoint[1] = v;
 671   Real rightGridPoint[2];
 672   rightGridPoint[0] = uright;
 673   rightGridPoint[1] = v;
 674
 675   Int i;
 676   Int index1, index2;
 677
 678   index1 = leftChain->findIndexFirstAboveEqualGen(v, leftChainStartIndex, leftChainEndIndex);
 679
 680
 681   index2 = rightChain->findIndexFirstAboveEqualGen(v, rightChainStartIndex, rightChainEndIndex);
 682
 683   if(index2>= leftChainStartIndex)  //index2 was found above
 684     index2 = rightChain->skipEqualityFromStart(v, index2, rightChainEndIndex);
 685
 686   if(index1<leftChainStartIndex && index2 <rightChainStartIndex) /*no point above v on left chain or right chain*/
 687     {
 688       /*the topVertex is the only vertex above v*/
 689       ret_leftCornerWhere = 1;
 690       ret_rightCornerWhere = 1;
 691     }
 692   else if(index1<leftChainStartIndex ) /*index2 >= rightChainStartIndex*/
 693     {
 694       ret_rightCornerWhere = 2; /*on right chain*/
 695       ret_rightCornerIndex = index2;
 696
 697       //find the minimum u on right top, either that, or top, or right[index2] is the left corner
 698       Real tempMin = rightChain->getVertex(index2)[0];
 699       Int tempI = index2;
 700       for(i=index2-1; i>=rightChainStartIndex; i--)
 701         if(rightChain->getVertex(i)[0] < tempMin)
 702           {
 703             tempMin = rightChain->getVertex(i)[0];
 704             tempI = i;
 705           }
 706       //chech whether (leftGridPoint, top) intersects rightchai,
 707       //if yes, use right corner as left corner
 708       //if not, use top or right[tempI] as left corner
 709       if(DBG_intersectChain(rightChain, rightChainStartIndex, rightChainEndIndex,
 710                         leftGridPoint, topVertex))
 711         {
 712           ret_leftCornerWhere = 2; //rightChain
 713           ret_leftCornerIndex = index2;
 714         }
 715       else if(topVertex[0] < tempMin)
 716         ret_leftCornerWhere = 1; /*topvertex*/
 717       else
 718         {
 719           ret_leftCornerWhere = 2; //right chain
 720           ret_leftCornerIndex = tempI;
 721         }
 722
 723     }
 724   else if(index2< rightChainStartIndex) /*index1>=leftChainStartIndex*/
 725     {
 726       ret_leftCornerWhere = 0; /*left chain*/
 727       ret_leftCornerIndex = index1;
 728
 729       //find the maximum u on the left top section. either that or topvertex, or left[index1]  is the right corner
 730       Real tempMax = leftChain->getVertex(index1)[0];
 731       Int tempI = index1;
 732
 733       for(i=index1-1; i>=leftChainStartIndex; i--){
 734
 735         if(leftChain->getVertex(i)[0] > tempMax)
 736           {
 737
 738             tempMax = leftChain->getVertex(i)[0];
 739             tempI = i;
 740           }
 741       }
 742       //check whether (rightGridPoint, top) intersects leftChain or not
 743       //if yes, we use leftCorner as the right corner
 744       //if not, we use either top or left[tempI] as the right corner
 745       if(DBG_intersectChain(leftChain, leftChainStartIndex,leftChainEndIndex,
 746                             rightGridPoint, topVertex))
 747          {
 748            ret_rightCornerWhere = 0; //left chan
 749            ret_rightCornerIndex = index1;
 750          }
 751       else if(topVertex[0] > tempMax)
 752         ret_rightCornerWhere = 1;//topVertex
 753       else
 754         {
 755           ret_rightCornerWhere = 0;//left chain
 756           ret_rightCornerIndex = tempI;
 757         }
 758     }
 759   else /*index1>=leftChainStartIndex and index2 >=rightChainStartIndex*/
 760     {
 761       if(leftChain->getVertex(index1)[1] <= rightChain->getVertex(index2)[1]) /*left point below right point*/
 762         {
 763           ret_leftCornerWhere = 0; /*on left chain*/
 764           ret_leftCornerIndex = index1;
 765
 766           Real tempMax;
 767           Int tempI;
 768
 769           tempI = index1;
 770           tempMax = leftChain->getVertex(index1)[0];
 771
 772           /*find the maximum u for all the points on the left below the right point index2*/
 773           for(i=index1-1; i>= leftChainStartIndex; i--)
 774             {
 775               if(leftChain->getVertex(i)[1] > rightChain->getVertex(index2)[1])
 776                 break;
 777
 778               if(leftChain->getVertex(i)[0]>tempMax)
 779                 {
 780                   tempI = i;
 781                   tempMax = leftChain->getVertex(i)[0];
 782                 }
 783             }
 784           //chek whether (rightChain(index2), rightGridPoint) intersects leftchian or not
 785           if(DBG_intersectChain(leftChain, leftChainStartIndex, leftChainEndIndex, rightGridPoint, rightChain->getVertex(index2)))
 786              {
 787                ret_rightCornerWhere = 0;
 788                ret_rightCornerIndex = index1;
 789              }
 790           else if(tempMax >= rightChain->getVertex(index2)[0] ||
 791              tempMax >= uright)
 792             {
 793               ret_rightCornerWhere = 0; /*on left Chain*/
 794               ret_rightCornerIndex = tempI;
 795             }
 796           else
 797             {
 798               ret_rightCornerWhere = 2; /*on right chain*/
 799               ret_rightCornerIndex = index2;
 800             }
 801         }
 802       else /*left above right*/
 803         {
 804           ret_rightCornerWhere = 2; /*on the right*/
 805           ret_rightCornerIndex = index2;
 806
 807           Real tempMin;
 808           Int tempI;
 809
 810           tempI = index2;
 811           tempMin = rightChain->getVertex(index2)[0];
 812
 813           /*find the minimum u for all the points on the right below the left poitn index1*/
 814           for(i=index2-1; i>= rightChainStartIndex; i--)
 815             {
 816               if( rightChain->getVertex(i)[1] > leftChain->getVertex(index1)[1])
 817                 break;
 818               if(rightChain->getVertex(i)[0] < tempMin)
 819                 {
 820                   tempI = i;
 821                   tempMin = rightChain->getVertex(i)[0];
 822                 }
 823             }
 824           //check whether (leftGRidPoint,left(index1)) interesect right chain
 825           if(DBG_intersectChain(rightChain, rightChainStartIndex, rightChainEndIndex,
 826                                 leftGridPoint, leftChain->getVertex(index1)))
 827             {
 828               ret_leftCornerWhere = 2; //right
 829               ret_leftCornerIndex = index2;
 830             }
 831           else if(tempMin <= leftChain->getVertex(index1)[0] ||
 832              tempMin <= uleft)
 833             {
 834               ret_leftCornerWhere = 2; /* on right chain*/
 835               ret_leftCornerIndex = tempI;
 836             }
 837           else
 838             {
 839               ret_leftCornerWhere = 0; /*on left chain*/
 840               ret_leftCornerIndex = index1;
 841             }
 842         }
 843     }
 844 #ifdef MYDEBUG
 845 printf("***********leave findUpCorners\n");
 846 #endif
 847 }
 848
 849 //return 1 if neck exists, 0 othewise
 850 Int findNeckF(vertexArray *leftChain, Int botLeftIndex,
 851               vertexArray *rightChain, Int botRightIndex,
 852               gridBoundaryChain* leftGridChain,
 853               gridBoundaryChain* rightGridChain,
 854               Int gridStartIndex,
 855               Int& neckLeft,
 856               Int& neckRight)
 857 {
 858 /*
 859 printf("enter findNeckF, botleft, botright=%i,%i,gstartindex=%i\n",botLeftIndex,botRightIndex,gridStartIndex);
 860 printf("leftChain is\n");
 861 leftChain->print();
 862 printf("rightChain is\n");
 863 rightChain->print();
 864 */
 865
 866   Int lowerGridIndex; //the grid below leftChain and rightChian vertices
 867   Int i;
 868   Int n_vlines = leftGridChain->get_nVlines();
 869   Real v;
 870   if(botLeftIndex >= leftChain->getNumElements() ||
 871      botRightIndex >= rightChain->getNumElements())
 872     return 0; //no neck exists
 873
 874   v=min(leftChain->getVertex(botLeftIndex)[1], rightChain->getVertex(botRightIndex)[1]);
 875
 876
 877
 878
 879   for(i=gridStartIndex; i<n_vlines; i++)
 880     if(leftGridChain->get_v_value(i) <= v &&
 881        leftGridChain->getUlineIndex(i)<= rightGridChain->getUlineIndex(i))
 882       break;
 883
 884   lowerGridIndex = i;
 885
 886   if(lowerGridIndex == n_vlines) //the two trm vertex are higher than all gridlines
 887     return 0;
 888   else
 889     {
 890       Int botLeft2, botRight2;
 891 /*
 892 printf("leftGridChain->get_v_)value=%f\n",leftGridChain->get_v_value(lowerGridIndex), botLeftIndex);
 893 printf("leftChain->get_vertex(0)=(%f,%f)\n", leftChain->getVertex(0)[0],leftChain->getVertex(0)[1]);
 894 printf("leftChain->get_vertex(1)=(%f,%f)\n", leftChain->getVertex(1)[0],leftChain->getVertex(1)[1]);
 895 printf("leftChain->get_vertex(2)=(%f,%f)\n", leftChain->getVertex(2)[0],leftChain->getVertex(2)[1]);
 896 */
 897       botLeft2 = leftChain->findIndexFirstAboveEqualGen(leftGridChain->get_v_value(lowerGridIndex), botLeftIndex, leftChain->getNumElements()-1) -1 ;
 898
 899 /*
 900 printf("botLeft2=%i\n", botLeft2);
 901 printf("leftChain->getNumElements=%i\n", leftChain->getNumElements());
 902 */
 903
 904       botRight2 = rightChain->findIndexFirstAboveEqualGen(leftGridChain->get_v_value(lowerGridIndex), botRightIndex, rightChain->getNumElements()-1) -1;
 905       if(botRight2 < botRightIndex) botRight2=botRightIndex;
 906
 907       if(botLeft2 < botLeftIndex) botLeft2 = botLeftIndex;
 908
 909       assert(botLeft2 >= botLeftIndex);
 910       assert(botRight2 >= botRightIndex);
 911       //find nectLeft so that it is th erightmost vertex on letChain
 912
 913       Int tempI = botLeftIndex;
 914       Real temp = leftChain->getVertex(tempI)[0];
 915       for(i=botLeftIndex+1; i<= botLeft2; i++)
 916         if(leftChain->getVertex(i)[0] > temp)
 917           {
 918             temp = leftChain->getVertex(i)[0];
 919             tempI = i;
 920           }
 921       neckLeft = tempI;
 922
 923       tempI = botRightIndex;
 924       temp = rightChain->getVertex(tempI)[0];
 925       for(i=botRightIndex+1; i<= botRight2; i++)
 926         if(rightChain->getVertex(i)[0] < temp)
 927           {
 928             temp = rightChain->getVertex(i)[0];
 929             tempI = i;
 930           }
 931       neckRight = tempI;
 932       return 1;
 933     }
 934 }
 935
 936
 937
 938 /*find i>=botLeftIndex,j>=botRightIndex so that
 939  *(leftChain[i], rightChain[j]) is a neck.
 940  */
 941 void findNeck(vertexArray *leftChain, Int botLeftIndex,
 942               vertexArray *rightChain, Int botRightIndex,
 943               Int& leftLastIndex, /*left point of the neck*/
 944               Int& rightLastIndex /*right point of the neck*/
 945               )
 946 {
 947   assert(botLeftIndex < leftChain->getNumElements() &&
 948      botRightIndex < rightChain->getNumElements());
 949
 950   /*now the neck exists for sure*/
 951
 952   if(leftChain->getVertex(botLeftIndex)[1] <= rightChain->getVertex(botRightIndex)[1]) //left below right
 953     {
 954
 955       leftLastIndex = botLeftIndex;
 956
 957       /*find i so that rightChain[i][1] >= leftchainbotverte[1], and i+1<
 958        */
 959       rightLastIndex=rightChain->findIndexAboveGen(leftChain->getVertex(botLeftIndex)[1], botRightIndex+1, rightChain->getNumElements()-1);
 960     }
 961   else  //left above right
 962     {
 963
 964       rightLastIndex = botRightIndex;
 965
 966       leftLastIndex = leftChain->findIndexAboveGen(rightChain->getVertex(botRightIndex)[1],
 967                                                   botLeftIndex+1,
 968                                                   leftChain->getNumElements()-1);
 969     }
 970 }
 971
 972
 973
 974 void findLeftGridIndices(directedLine* topEdge, Int firstGridIndex, Int lastGridIndex, gridWrap* grid,  Int* ret_indices, Int* ret_innerIndices)
 975 {
 976
 977   Int i,k,isHoriz = 0;
 978   Int n_ulines = grid->get_n_ulines();
 979   Real uMin = grid->get_u_min();
 980   Real uMax = grid->get_u_max();
 981   /*
 982   Real vMin = grid->get_v_min();
 983   Real vMax = grid->get_v_max();
 984   */
 985   Real slop = 0.0, uinterc;
 986
 987 #ifdef SHORTEN_GRID_LINE
 988   //uintercBuf stores all the interction u value for each grid line
 989   //notice that lastGridIndex<= firstGridIndex
 990   Real *uintercBuf = (Real *) malloc (sizeof(Real) * (firstGridIndex-lastGridIndex+1));
 991   assert(uintercBuf);
 992 #endif
 993
 994   /*initialization to make vtail bigger than grid->...*/
 995   directedLine* dLine = topEdge;
 996   Real vtail = grid->get_v_value(firstGridIndex) + 1.0;
 997   Real tempMaxU = grid->get_u_min();
 998
 999
1000   /*for each grid line*/
1001   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1002     {
1003
1004       Real grid_v_value = grid->get_v_value(i);
1005
1006       /*check whether this grid line is below the current trim edge.*/
1007       if(vtail > grid_v_value)
1008         {
1009           /*since the grid line is below the trim edge, we
1010            *find the trim edge which will contain the trim line
1011            */
1012           while( (vtail=dLine->tail()[1]) > grid_v_value){
1013
1014             tempMaxU = max(tempMaxU, dLine->tail()[0]);
1015             dLine = dLine -> getNext();
1016           }
1017
1018           if( fabs(dLine->head()[1] - vtail) < ZERO)
1019             isHoriz = 1;
1020           else
1021             {
1022               isHoriz = 0;
1023               slop = (dLine->head()[0] - dLine->tail()[0]) / (dLine->head()[1]-vtail);
1024             }
1025         }
1026
1027       if(isHoriz)
1028         {
1029           uinterc = max(dLine->head()[0], dLine->tail()[0]);
1030         }
1031       else
1032         {
1033           uinterc = slop * (grid_v_value - vtail) + dLine->tail()[0];
1034         }
1035
1036       tempMaxU = max(tempMaxU, uinterc);
1037
1038       if(uinterc < uMin && uinterc >= uMin - ZERO)
1039         uinterc = uMin;
1040       if(uinterc > uMax && uinterc <= uMax + ZERO)
1041         uinterc = uMax;
1042
1043 #ifdef SHORTEN_GRID_LINE
1044       uintercBuf[k] = uinterc;
1045 #endif
1046
1047       assert(uinterc >= uMin && uinterc <= uMax);
1048        if(uinterc == uMax)
1049          ret_indices[k] = n_ulines-1;
1050        else
1051          ret_indices[k] = (Int)(((uinterc-uMin)/(uMax - uMin)) * (n_ulines-1)) + 1;
1052       if(ret_indices[k] >= n_ulines)
1053         ret_indices[k] = n_ulines-1;
1054
1055
1056       ret_innerIndices[k] = (Int)(((tempMaxU-uMin)/(uMax - uMin)) * (n_ulines-1)) + 1;
1057
1058       /*reinitialize tempMaxU for next grdiLine*/
1059       tempMaxU = uinterc;
1060     }
1061 #ifdef SHORTEN_GRID_LINE
1062   //for each grid line, compare the left grid point with the
1063   //intersection point. If the two points are too close, then
1064   //we should move the grid point one grid to the right
1065   //and accordingly we should update the inner index.
1066   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1067     {
1068       //check gridLine i
1069       //check ret_indices[k]
1070       Real a = grid->get_u_value(ret_indices[k]-1);
1071       Real b = grid->get_u_value(ret_indices[k]);
1072       assert(uintercBuf[k] >= a && uintercBuf < b);
1073       if( (b-uintercBuf[k]) <= 0.2 * (b-a)) //interc is very close to b
1074         {
1075           ret_indices[k]++;
1076         }
1077
1078       //check ret_innerIndices[k]
1079       if(k>0)
1080         {
1081           if(ret_innerIndices[k] < ret_indices[k-1])
1082             ret_innerIndices[k] = ret_indices[k-1];
1083           if(ret_innerIndices[k] < ret_indices[k])
1084             ret_innerIndices[k] = ret_indices[k];
1085         }
1086     }
1087   //clean up
1088   free(uintercBuf);
1089 #endif
1090 }
1091
1092 void findRightGridIndices(directedLine* topEdge, Int firstGridIndex, Int lastGridIndex, gridWrap* grid,  Int* ret_indices, Int* ret_innerIndices)
1093 {
1094
1095   Int i,k;
1096   Int n_ulines = grid->get_n_ulines();
1097   Real uMin = grid->get_u_min();
1098   Real uMax = grid->get_u_max();
1099   /*
1100   Real vMin = grid->get_v_min();
1101   Real vMax = grid->get_v_max();
1102   */
1103   Real slop = 0.0, uinterc;
1104
1105 #ifdef SHORTEN_GRID_LINE
1106   //uintercBuf stores all the interction u value for each grid line
1107   //notice that firstGridIndex >= lastGridIndex
1108   Real *uintercBuf = (Real *) malloc (sizeof(Real) * (firstGridIndex-lastGridIndex+1));
1109   assert(uintercBuf);
1110 #endif
1111
1112   /*initialization to make vhead bigger than grid->v_value...*/
1113   directedLine* dLine = topEdge->getPrev();
1114   Real vhead = dLine->tail()[1];
1115   Real tempMinU = grid->get_u_max();
1116
1117   /*for each grid line*/
1118   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1119     {
1120
1121       Real grid_v_value = grid->get_v_value(i);
1122
1123
1124       /*check whether this grid line is below the current trim edge.*/
1125       if(vhead >= grid_v_value)
1126         {
1127           /*since the grid line is below the tail of the trim edge, we
1128            *find the trim edge which will contain the trim line
1129            */
1130           while( (vhead=dLine->head()[1]) > grid_v_value){
1131             tempMinU = min(tempMinU, dLine->head()[0]);
1132             dLine = dLine -> getPrev();
1133           }
1134
1135           /*skip the equality in the case of degenerat case: horizontal */
1136           while(dLine->head()[1] == grid_v_value)
1137             dLine = dLine->getPrev();
1138
1139           assert( dLine->tail()[1] != dLine->head()[1]);
1140           slop = (dLine->tail()[0] - dLine->head()[0]) / (dLine->tail()[1]-dLine->head()[1]);
1141           /*
1142            if(dLine->tail()[1] == vhead)
1143              isHoriz = 1;
1144              else
1145             {
1146               isHoriz = 0;
1147               slop = (dLine->tail()[0] - dLine->head()[0]) / (dLine->tail()[1]-vhead);
1148             }
1149             */
1150         }
1151         uinterc = slop * (grid_v_value - dLine->head()[1]) + dLine->head()[0];
1152
1153       //in case unterc is outside of the grid due to floating point
1154       if(uinterc < uMin)
1155         uinterc = uMin;
1156       else if(uinterc > uMax)
1157         uinterc = uMax;
1158
1159 #ifdef SHORTEN_GRID_LINE
1160       uintercBuf[k] = uinterc;
1161 #endif
1162
1163       tempMinU = min(tempMinU, uinterc);
1164
1165       assert(uinterc >= uMin && uinterc <= uMax);
1166
1167       if(uinterc == uMin)
1168         ret_indices[k] = 0;
1169       else
1170         ret_indices[k] = (int)ceil((((uinterc-uMin)/(uMax - uMin)) * (n_ulines-1))) -1;
1171 /*
1172 if(ret_indices[k] >= grid->get_n_ulines())
1173   {
1174   printf("ERROR3\n");
1175   exit(0);
1176 }
1177 if(ret_indices[k] < 0)
1178   {
1179   printf("ERROR4\n");
1180   exit(0);
1181 }
1182 */
1183       ret_innerIndices[k] = (int)ceil ((((tempMinU-uMin)/(uMax - uMin)) * (n_ulines-1))) -1;
1184
1185       tempMinU = uinterc;
1186     }
1187 #ifdef SHORTEN_GRID_LINE
1188   //for each grid line, compare the left grid point with the
1189   //intersection point. If the two points are too close, then
1190   //we should move the grid point one grid to the right
1191   //and accordingly we should update the inner index.
1192   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1193     {
1194       //check gridLine i
1195       //check ret_indices[k]
1196       Real a = grid->get_u_value(ret_indices[k]);
1197       Real b = grid->get_u_value(ret_indices[k]+1);
1198       assert(uintercBuf[k] > a && uintercBuf <= b);
1199       if( (uintercBuf[k]-a) <= 0.2 * (b-a)) //interc is very close to a
1200         {
1201           ret_indices[k]--;
1202         }
1203
1204       //check ret_innerIndices[k]
1205       if(k>0)
1206         {
1207           if(ret_innerIndices[k] > ret_indices[k-1])
1208             ret_innerIndices[k] = ret_indices[k-1];
1209           if(ret_innerIndices[k] > ret_indices[k])
1210             ret_innerIndices[k] = ret_indices[k];
1211         }
1212     }
1213   //clean up
1214   free(uintercBuf);
1215 #endif
1216 }
1217
1218
1219 void sampleMonoPoly(directedLine* polygon, gridWrap* grid, Int ulinear, Int vlinear, primStream* pStream, rectBlockArray* rbArray)
1220 {
1221 /*
1222 {
1223 grid->print();
1224 polygon->writeAllPolygons("zloutputFile");
1225 exit(0);
1226 }
1227 */
1228
1229 if(grid->get_n_ulines() == 2 ||
1230    grid->get_n_vlines() == 2)
1231 {
1232   if(ulinear && grid->get_n_ulines() == 2)
1233     {
1234       monoTriangulationFun(polygon, compV2InY, pStream);
1235       return;
1236     }
1237   else if(DBG_isConvex(polygon) && polygon->numEdges() >=4)
1238      {
1239        triangulateConvexPoly(polygon, ulinear, vlinear, pStream);
1240        return;
1241      }
1242   else if(vlinear || DBG_is_U_direction(polygon))
1243     {
1244       Int n_cusps;//num interior cusps
1245       Int n_edges = polygon->numEdges();
1246       directedLine** cusps = (directedLine**) malloc(sizeof(directedLine*) * n_edges);
1247       assert(cusps);
1248       findInteriorCuspsX(polygon, n_cusps, cusps);
1249
1250       if(n_cusps == 0) //u_monotone
1251         {
1252
1253           monoTriangulationFun(polygon, compV2InX, pStream);
1254
1255           free(cusps);
1256           return;
1257         }
1258       else if(n_cusps == 1) //one interior cusp
1259         {
1260
1261           directedLine* new_polygon = polygonConvert(cusps[0]);
1262
1263           directedLine* other = findDiagonal_singleCuspX( new_polygon);
1264
1265
1266
1267           //<other> should NOT be null unless there are self-intersecting
1268           //trim curves. In that case, we don't want to core dump, instead,
1269           //we triangulate anyway, and print out error message.
1270           if(other == NULL)
1271             {
1272               monoTriangulationFun(polygon, compV2InX, pStream);
1273               free(cusps);
1274               return;
1275             }
1276
1277           directedLine* ret_p1;
1278           directedLine* ret_p2;
1279
1280           new_polygon->connectDiagonal_2slines(new_polygon, other,
1281                                            &ret_p1,
1282                                            &ret_p2,
1283                                            new_polygon);
1284
1285           monoTriangulationFun(ret_p1, compV2InX, pStream);
1286           monoTriangulationFun(ret_p2, compV2InX, pStream);
1287
1288           ret_p1->deleteSinglePolygonWithSline();
1289           ret_p2->deleteSinglePolygonWithSline();
1290
1291           free(cusps);
1292           return;
1293          }
1294      free(cusps);
1295      }
1296 }
1297
1298   /*find the top and bottom of the polygon. It is supposed to be
1299    *a V-monotone polygon
1300    */
1301
1302   directedLine* tempV;
1303   directedLine* topV;
1304   directedLine* botV;
1305   topV = botV = polygon;
1306
1307   for(tempV = polygon->getNext(); tempV != polygon; tempV = tempV->getNext())
1308     {
1309       if(compV2InY(topV->head(), tempV->head())<0) {
1310
1311         topV = tempV;
1312       }
1313       if(compV2InY(botV->head(), tempV->head())>0) {
1314
1315         botV = tempV;
1316       }
1317     }
1318
1319   /*find the first(top) and the last (bottom) grid line which intersect the
1320    *this polygon
1321    */
1322   Int firstGridIndex; /*the index in the grid*/
1323   Int lastGridIndex;
1324   firstGridIndex = (Int) ((topV->head()[1] - grid->get_v_min()) / (grid->get_v_max() - grid->get_v_min()) * (grid->get_n_vlines()-1));
1325   lastGridIndex  = (Int) ((botV->head()[1] - grid->get_v_min()) / (grid->get_v_max() - grid->get_v_min()) * (grid->get_n_vlines()-1)) + 1;
1326
1327
1328   /*find the interval inside  the polygon for each gridline*/
1329   Int *leftGridIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1330   assert(leftGridIndices);
1331   Int *rightGridIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1332   assert(rightGridIndices);
1333   Int *leftGridInnerIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1334   assert(leftGridInnerIndices);
1335   Int *rightGridInnerIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1336   assert(rightGridInnerIndices);
1337
1338   findLeftGridIndices(topV, firstGridIndex, lastGridIndex, grid,  leftGridIndices, leftGridInnerIndices);
1339
1340   findRightGridIndices(topV, firstGridIndex, lastGridIndex, grid,  rightGridIndices, rightGridInnerIndices);
1341
1342   gridBoundaryChain leftGridChain(grid, firstGridIndex, firstGridIndex-lastGridIndex+1, leftGridIndices, leftGridInnerIndices);
1343
1344   gridBoundaryChain rightGridChain(grid, firstGridIndex, firstGridIndex-lastGridIndex+1, rightGridIndices, rightGridInnerIndices);
1345
1346
1347
1348 //  leftGridChain.draw();
1349 //  leftGridChain.drawInner();
1350 //  rightGridChain.draw();
1351 //  rightGridChain.drawInner();
1352   /*(1) determine the grid boundaries (left and right).
1353    *(2) process polygon  into two monotone chaines: use vertexArray
1354    *(3) call sampleMonoPolyRec
1355    */
1356
1357   /*copy the two chains into vertexArray datastructure*/
1358   Int i;
1359   vertexArray leftChain(20); /*this is a dynamic array*/
1360   for(i=1; i<=topV->get_npoints()-2; i++) { /*the first vertex is the top vertex which doesn't belong to inc_chain*/
1361     leftChain.appendVertex(topV->getVertex(i));
1362   }
1363   for(tempV = topV->getNext(); tempV != botV; tempV = tempV->getNext())
1364     {
1365       for(i=0; i<=tempV->get_npoints()-2; i++){
1366         leftChain.appendVertex(tempV->getVertex(i));
1367       }
1368     }
1369
1370   vertexArray rightChain(20);
1371   for(tempV = topV->getPrev(); tempV != botV; tempV = tempV->getPrev())
1372     {
1373       for(i=tempV->get_npoints()-2; i>=0; i--){
1374         rightChain.appendVertex(tempV->getVertex(i));
1375       }
1376     }
1377   for(i=botV->get_npoints()-2; i>=1; i--){
1378     rightChain.appendVertex(tempV->getVertex(i));
1379   }
1380
1381   sampleMonoPolyRec(topV->head(),
1382                     botV->head(),
1383                     &leftChain,
1384                     0,
1385                     &rightChain,
1386                     0,
1387                     &leftGridChain,
1388                     &rightGridChain,
1389                     0,
1390                     pStream,
1391                     rbArray);
1392
1393
1394   /*cleanup space*/
1395   free(leftGridIndices);
1396   free(rightGridIndices);
1397   free(leftGridInnerIndices);
1398   free(rightGridInnerIndices);
1399 }
1400
1401 void sampleMonoPolyRec(
1402                        Real* topVertex,
1403                        Real* botVertex,
1404                        vertexArray* leftChain,
1405                        Int leftStartIndex,
1406                        vertexArray* rightChain,
1407                        Int rightStartIndex,
1408                        gridBoundaryChain* leftGridChain,
1409                        gridBoundaryChain* rightGridChain,
1410                        Int gridStartIndex,
1411                        primStream* pStream,
1412                        rectBlockArray* rbArray)
1413 {
1414
1415   /*find the first connected component, and the four corners.
1416    */
1417   Int index1, index2; /*the first and last grid line of the first connected component*/
1418
1419   if(topVertex[1] <= botVertex[1])
1420     return;
1421
1422   /*find i so that the grid line is below the top vertex*/
1423   Int i=gridStartIndex;
1424   while (i < leftGridChain->get_nVlines())
1425     {
1426       if(leftGridChain->get_v_value(i) < topVertex[1])
1427         break;
1428       i++;
1429     }
1430
1431   /*find the first connected component starting with i*/
1432   /*find index1 so that left_uline_index <= right_uline_index, that is, this
1433    *grid line contains at least one inner grid point
1434    */
1435   index1=i;
1436   int num_skipped_grid_lines=0;
1437   while(index1 < leftGridChain->get_nVlines())
1438     {
1439       if(leftGridChain->getUlineIndex(index1) <= rightGridChain->getUlineIndex(index1))
1440         break;
1441       num_skipped_grid_lines++;
1442       index1++;
1443     }
1444
1445
1446
1447   if(index1 >= leftGridChain->get_nVlines()) /*no grid line exists which has inner point*/
1448     {
1449       /*stop recursion, ...*/
1450       /*monotone triangulate it...*/
1451 //      printf("no grid line exists\n");
1452 /*
1453       monoTriangulationRecOpt(topVertex, botVertex, leftChain, leftStartIndex,
1454                            rightChain, rightStartIndex, pStream);
1455 */
1456
1457 if(num_skipped_grid_lines <2)
1458   {
1459     monoTriangulationRecGenOpt(topVertex, botVertex, leftChain, leftStartIndex,
1460                                leftChain->getNumElements()-1,
1461                                rightChain, rightStartIndex,
1462                                rightChain->getNumElements()-1,
1463                                pStream);
1464   }
1465 else
1466   {
1467     //the optimum way to triangulate is top-down since this polygon
1468     //is narrow-long.
1469     monoTriangulationRec(topVertex, botVertex, leftChain, leftStartIndex,
1470                          rightChain, rightStartIndex, pStream);
1471   }
1472
1473 /*
1474       monoTriangulationRec(topVertex, botVertex, leftChain, leftStartIndex,
1475                            rightChain, rightStartIndex, pStream);
1476 */
1477
1478 /*      monoTriangulationRecGenTBOpt(topVertex, botVertex,
1479                                    leftChain, leftStartIndex, leftChain->getNumElements()-1,
1480                                    rightChain, rightStartIndex, rightChain->getNumElements()-1,
1481                                    pStream);*/
1482
1483
1484
1485     }
1486   else
1487     {
1488
1489       /*find index2 so that left_inner_index <= right_inner_index holds until index2*/
1490       index2=index1+1;
1491       if(index2 < leftGridChain->get_nVlines())
1492         while(leftGridChain->getInnerIndex(index2) <= rightGridChain->getInnerIndex(index2))
1493           {
1494             index2++;
1495             if(index2 >= leftGridChain->get_nVlines())
1496               break;
1497           }
1498
1499       index2--;
1500
1501
1502
1503       /*the neck*/
1504       Int neckLeftIndex;
1505       Int neckRightIndex;
1506
1507       /*the four corners*/
1508       Int up_leftCornerWhere;
1509       Int up_leftCornerIndex;
1510       Int up_rightCornerWhere;
1511       Int up_rightCornerIndex;
1512       Int down_leftCornerWhere;
1513       Int down_leftCornerIndex;
1514       Int down_rightCornerWhere;
1515       Int down_rightCornerIndex;
1516
1517       Real* tempBotVertex; /*the bottom vertex for this component*/
1518       Real* nextTopVertex=NULL; /*for the recursion*/
1519       Int nextLeftStartIndex=0;
1520       Int nextRightStartIndex=0;
1521
1522       /*find the points below the grid line index2 on both chains*/
1523       Int botLeftIndex = leftChain->findIndexStrictBelowGen(
1524                                                       leftGridChain->get_v_value(index2),
1525                                                       leftStartIndex,
1526                                                       leftChain->getNumElements()-1);
1527       Int botRightIndex = rightChain->findIndexStrictBelowGen(
1528                                                         rightGridChain->get_v_value(index2),
1529                                                         rightStartIndex,
1530                                                         rightChain->getNumElements()-1);
1531       /*if either botLeftIndex>= numelements,
1532        *        or botRightIndex >= numelemnet,
1533        *then there is no neck exists. the bottom vertex is botVertex,
1534        */
1535       if(! findNeckF(leftChain, botLeftIndex, rightChain, botRightIndex,
1536                    leftGridChain, rightGridChain, index2, neckLeftIndex, neckRightIndex))
1537          /*
1538       if(botLeftIndex == leftChain->getNumElements() ||
1539          botRightIndex == rightChain->getNumElements())
1540          */
1541         {
1542 #ifdef MYDEBUG
1543           printf("neck NOT exists, botRightIndex=%i\n", botRightIndex);
1544 #endif
1545
1546           tempBotVertex =  botVertex;
1547           nextTopVertex = botVertex;
1548           botLeftIndex = leftChain->getNumElements()-1;
1549           botRightIndex = rightChain->getNumElements()-1;
1550         }
1551       else /*neck exists*/
1552         {
1553 #ifdef MYDEBUG
1554           printf("neck exists\n");
1555 #endif
1556
1557           /*
1558           findNeck(leftChain, botLeftIndex,
1559                    rightChain, botRightIndex,
1560                    neckLeftIndex,
1561                    neckRightIndex);
1562                    */
1563 #ifdef MYDEBUG
1564 printf("neck is found, neckLeftIndex=%i, neckRightIndex=%i\n", neckLeftIndex, neckRightIndex);
1565 glBegin(GL_LINES);
1566 glVertex2fv(leftChain->getVertex(neckLeftIndex));
1567 glVertex2fv(rightChain->getVertex(neckRightIndex));
1568 glEnd();
1569 #endif
1570
1571           if(leftChain->getVertex(neckLeftIndex)[1] <= rightChain->getVertex(neckRightIndex)[1])
1572             {
1573               tempBotVertex = leftChain->getVertex(neckLeftIndex);
1574               botLeftIndex = neckLeftIndex-1;
1575               botRightIndex = neckRightIndex;
1576               nextTopVertex = rightChain->getVertex(neckRightIndex);
1577               nextLeftStartIndex = neckLeftIndex;
1578               nextRightStartIndex = neckRightIndex+1;
1579             }
1580           else
1581             {
1582               tempBotVertex = rightChain->getVertex(neckRightIndex);
1583               botLeftIndex = neckLeftIndex;
1584               botRightIndex = neckRightIndex-1;
1585               nextTopVertex = leftChain->getVertex(neckLeftIndex);
1586               nextLeftStartIndex = neckLeftIndex+1;
1587               nextRightStartIndex = neckRightIndex;
1588             }
1589         }
1590
1591       findUpCorners(topVertex,
1592                     leftChain,
1593                     leftStartIndex, botLeftIndex,
1594                     rightChain,
1595                     rightStartIndex, botRightIndex,
1596                     leftGridChain->get_v_value(index1),
1597                     leftGridChain->get_u_value(index1),
1598                     rightGridChain->get_u_value(index1),
1599                     up_leftCornerWhere,
1600                     up_leftCornerIndex,
1601                     up_rightCornerWhere,
1602                     up_rightCornerIndex);
1603
1604       findDownCorners(tempBotVertex,
1605                       leftChain,
1606                       leftStartIndex, botLeftIndex,
1607                       rightChain,
1608                       rightStartIndex, botRightIndex,
1609                       leftGridChain->get_v_value(index2),
1610                       leftGridChain->get_u_value(index2),
1611                       rightGridChain->get_u_value(index2),
1612                       down_leftCornerWhere,
1613                       down_leftCornerIndex,
1614                       down_rightCornerWhere,
1615                       down_rightCornerIndex);
1616 #ifdef MYDEBUG
1617       printf("find corners done, down_leftwhere=%i, down_righwhere=%i,\n",down_leftCornerWhere, down_rightCornerWhere );
1618       printf("find corners done, up_leftwhere=%i, up_righwhere=%i,\n",up_leftCornerWhere, up_rightCornerWhere );
1619       printf("find corners done, up_leftindex=%i, up_righindex=%i,\n",up_leftCornerIndex, up_rightCornerIndex );
1620       printf("find corners done, down_leftindex=%i, down_righindex=%i,\n",down_leftCornerIndex, down_rightCornerIndex );
1621 #endif
1622
1623 /*
1624       drawCorners(topVertex,
1625                   tempBotVertex,
1626                   leftChain,
1627                   rightChain,
1628                   leftGridChain,
1629                   rightGridChain,
1630                   index1,
1631                   index2,
1632                   up_leftCornerWhere,
1633                   up_leftCornerIndex,
1634                   up_rightCornerWhere,
1635                   up_rightCornerIndex,
1636                   down_leftCornerWhere,
1637                   down_leftCornerIndex,
1638                   down_rightCornerWhere,
1639                   down_rightCornerIndex);
1640 */
1641
1642
1643       sampleConnectedComp(topVertex, tempBotVertex,
1644                           leftChain,
1645                           leftStartIndex, botLeftIndex,
1646                           rightChain,
1647                           rightStartIndex, botRightIndex,
1648                           leftGridChain,
1649                           rightGridChain,
1650                           index1, index2,
1651                           up_leftCornerWhere,
1652                           up_leftCornerIndex,
1653                           up_rightCornerWhere,
1654                           up_rightCornerIndex,
1655                           down_leftCornerWhere,
1656                           down_leftCornerIndex,
1657                           down_rightCornerWhere,
1658                           down_rightCornerIndex,
1659                           pStream,
1660                           rbArray
1661                           );
1662
1663       /*recursion*/
1664
1665       sampleMonoPolyRec(
1666                         nextTopVertex,
1667                         botVertex,
1668                         leftChain,
1669                         nextLeftStartIndex,
1670                         rightChain,
1671                         nextRightStartIndex,
1672                         leftGridChain,
1673                         rightGridChain,
1674                         index2+1,
1675                         pStream, rbArray);
1676
1677
1678     }
1679
1680 }
1681
1682 void sampleLeftStrip(vertexArray* leftChain,
1683                      Int topLeftIndex,
1684                      Int botLeftIndex,
1685                      gridBoundaryChain* leftGridChain,
1686                      Int leftGridChainStartIndex,
1687                      Int leftGridChainEndIndex,
1688                      primStream* pStream
1689                      )
1690 {
1691   assert(leftChain->getVertex(topLeftIndex)[1] > leftGridChain->get_v_value(leftGridChainStartIndex));
1692   assert(leftChain->getVertex(topLeftIndex+1)[1] <= leftGridChain->get_v_value(leftGridChainStartIndex));
1693   assert(leftChain->getVertex(botLeftIndex)[1] <= leftGridChain->get_v_value(leftGridChainEndIndex));
1694   assert(leftChain->getVertex(botLeftIndex-1)[1] > leftGridChain->get_v_value(leftGridChainEndIndex));
1695
1696   /*
1697    *(1)find the last grid line which doesn'; pass below
1698    * this first edge, sample this region: one trim edge and
1699    * possily multiple grid lines.
1700    */
1701   Real *upperVert, *lowerVert; /*the end points of the first trim edge*/
1702   upperVert = leftChain->getVertex(topLeftIndex);
1703   lowerVert = leftChain->getVertex(topLeftIndex+1);
1704
1705   Int index = leftGridChainStartIndex;
1706   while(leftGridChain->get_v_value(index) >= lowerVert[1]){
1707     index++;
1708     if(index > leftGridChainEndIndex)
1709       break;
1710   }
1711   index--;
1712
1713   sampleLeftSingleTrimEdgeRegion(upperVert, lowerVert,
1714                                  leftGridChain,
1715                                  leftGridChainStartIndex,
1716                                  index,
1717                                  pStream);
1718   sampleLeftStripRec(leftChain, topLeftIndex+1, botLeftIndex,
1719                      leftGridChain, index, leftGridChainEndIndex,
1720                      pStream);
1721
1722 }
1723
1724 void sampleLeftStripRec(vertexArray* leftChain,
1725                      Int topLeftIndex,
1726                      Int botLeftIndex,
1727                      gridBoundaryChain* leftGridChain,
1728                      Int leftGridChainStartIndex,
1729                      Int leftGridChainEndIndex,
1730                         primStream* pStream
1731                      )
1732 {
1733   /*now top left trim vertex is below the top grid line.
1734    */
1735   /*stop condition: if topLeftIndex >= botLeftIndex, then stop.
1736    */
1737   if(topLeftIndex >= botLeftIndex)
1738     return;
1739
1740   /*find the last trim vertex which is above the second top grid line:
1741    * index1.
1742    *and sampleLeftOneGridStep(leftchain, topLeftIndex, index1, leftGridChain,
1743    * leftGridChainStartIndex).
1744    * index1 could be equal to topLeftIndex.
1745    */
1746   Real secondGridChainV = leftGridChain->get_v_value(leftGridChainStartIndex+1);
1747   assert(leftGridChainStartIndex < leftGridChainEndIndex);
1748   Int index1 = topLeftIndex;
1749   while(leftChain->getVertex(index1)[1] > secondGridChainV)
1750     index1++;
1751   index1--;
1752
1753   sampleLeftOneGridStep(leftChain, topLeftIndex, index1, leftGridChain, leftGridChainStartIndex, pStream);
1754
1755
1756   /*
1757    * Let the next trim vertex be nextTrimVertIndex (which should be
1758    *  below the second grid line).
1759    * Find the last grid line index2 which is above nextTrimVert.
1760    * sampleLeftSingleTrimEdgeRegion(uppervert[2], lowervert[2],
1761    *                      leftGridChain, leftGridChainStartIndex+1, index2).
1762    */
1763   Real *uppervert, *lowervert;
1764   uppervert = leftChain->getVertex(index1);
1765   lowervert = leftChain->getVertex(index1+1);
1766   Int index2 = leftGridChainStartIndex+1;
1767
1768   while(leftGridChain->get_v_value(index2) >= lowervert[1])
1769     {
1770       index2++;
1771       if(index2 > leftGridChainEndIndex)
1772         break;
1773     }
1774   index2--;
1775   sampleLeftSingleTrimEdgeRegion(uppervert, lowervert, leftGridChain, leftGridChainStartIndex+1, index2,  pStream);
1776
1777    /* sampleLeftStripRec(leftChain,
1778                         nextTrimVertIndex,
1779                         botLeftIndex,
1780                         leftGridChain,
1781                         index2,
1782                         leftGridChainEndIndex
1783                         )
1784    *
1785    */
1786   sampleLeftStripRec(leftChain, index1+1, botLeftIndex, leftGridChain, index2, leftGridChainEndIndex, pStream);
1787
1788 }
1789
1790
1791 /***************begin RecF***********************/
1792 /* the gridlines from leftGridChainStartIndex to
1793  * leftGridChainEndIndex are assumed to form a
1794  * connected component.
1795  * the trim vertex of topLeftIndex is assumed to
1796  * be below the first gridline, and the tim vertex
1797  * of botLeftIndex is assumed to be above the last
1798  * grid line.
1799  * If botLeftIndex < topLeftIndex, then no connected componeent exists, and this funcion returns without
1800  * outputing any triangles.
1801  * Otherwise botLeftIndex >= topLeftIndex, there is at least one triangle to output.
1802  */
1803 void sampleLeftStripRecF(vertexArray* leftChain,
1804                      Int topLeftIndex,
1805                      Int botLeftIndex,
1806                      gridBoundaryChain* leftGridChain,
1807                      Int leftGridChainStartIndex,
1808                      Int leftGridChainEndIndex,
1809                         primStream* pStream
1810                      )
1811 {
1812   /*now top left trim vertex is below the top grid line.
1813    */
1814   /*stop condition: if topLeftIndex > botLeftIndex, then stop.
1815    */
1816   if(topLeftIndex > botLeftIndex)
1817     return;
1818
1819   /*if there is only one grid Line, return.*/
1820
1821   if(leftGridChainStartIndex>=leftGridChainEndIndex)
1822     return;
1823
1824
1825   assert(leftChain->getVertex(topLeftIndex)[1] <= leftGridChain->get_v_value(leftGridChainStartIndex) &&
1826          leftChain->getVertex(botLeftIndex)[1] >= leftGridChain->get_v_value(leftGridChainEndIndex));
1827
1828   /*firs find the first trim vertex which is below or equal to the second top grid line:
1829    * index1.
1830    */
1831   Real secondGridChainV = leftGridChain->get_v_value(leftGridChainStartIndex+1);
1832
1833
1834   Int index1 = topLeftIndex;
1835
1836   while(leftChain->getVertex(index1)[1] > secondGridChainV){
1837     index1++;
1838     if(index1>botLeftIndex)
1839       break;
1840   }
1841
1842   /*now leftChain->getVertex(index-1)[1] > secondGridChainV and
1843    *    leftChain->getVertex(index)[1] <= secondGridChainV
1844    *If equality holds, then we should include the vertex index1, otherwise we include only index1-1, to
1845    *perform sampleOneGridStep.
1846    */
1847   if(index1>botLeftIndex)
1848     index1--;
1849   else if(leftChain->getVertex(index1)[1] < secondGridChainV)
1850     index1--;
1851
1852   /*now we have leftChain->getVertex(index1)[1] >= secondGridChainV, and
1853    *  leftChain->getVertex(index1+1)[1] <= secondGridChainV
1854    */
1855
1856
1857   sampleLeftOneGridStep(leftChain, topLeftIndex, index1, leftGridChain, leftGridChainStartIndex, pStream);
1858
1859
1860   /*if leftChain->getVertex(index1)[1] == secondGridChainV, then we can recursively do the rest.
1861    */
1862   if(leftChain->getVertex(index1)[1] == secondGridChainV)
1863     {
1864
1865       sampleLeftStripRecF(leftChain, index1, botLeftIndex,leftGridChain, leftGridChainStartIndex+1, leftGridChainEndIndex, pStream);
1866     }
1867   else if(index1 < botLeftIndex)
1868     {
1869
1870       /* Otherwise, we have leftChain->getVertex(index1)[1] > secondGridChainV,
1871        * let the next trim vertex be nextTrimVertIndex (which should be  strictly
1872        *  below the second grid line).
1873        * Find the last grid line index2 which is above nextTrimVert.
1874        * sampleLeftSingleTrimEdgeRegion(uppervert[2], lowervert[2],
1875        *                      leftGridChain, leftGridChainStartIndex+1, index2).
1876        */
1877       Real *uppervert, *lowervert;
1878       uppervert = leftChain->getVertex(index1);
1879       lowervert = leftChain->getVertex(index1+1); //okay since index1<botLeftIndex
1880       Int index2 = leftGridChainStartIndex+1;
1881
1882
1883       while(leftGridChain->get_v_value(index2) >= lowervert[1])
1884         {
1885           index2++;
1886           if(index2 > leftGridChainEndIndex)
1887             break;
1888         }
1889       index2--;
1890
1891
1892       sampleLeftSingleTrimEdgeRegion(uppervert, lowervert, leftGridChain, leftGridChainStartIndex+1, index2,  pStream);
1893
1894       /*recursion*/
1895
1896       sampleLeftStripRecF(leftChain, index1+1, botLeftIndex, leftGridChain, index2, leftGridChainEndIndex, pStream);
1897     }
1898
1899 }
1900
1901 /***************End RecF***********************/
1902
1903 /*sample the left area in between one trim edge and multiple grid lines.
1904  * all the grid lines should be in between the two end poins of the
1905  *trim edge.
1906  */
1907 void sampleLeftSingleTrimEdgeRegion(Real upperVert[2], Real lowerVert[2],
1908                                     gridBoundaryChain* gridChain,
1909                                     Int beginIndex,
1910                                     Int endIndex,
1911                                     primStream* pStream)
1912 {
1913   Int i,j,k;
1914
1915   vertexArray vArray(endIndex-beginIndex+1);
1916   vArray.appendVertex(gridChain->get_vertex(beginIndex));
1917
1918   for(k=1, i=beginIndex+1; i<=endIndex; i++, k++)
1919     {
1920       vArray.appendVertex(gridChain->get_vertex(i));
1921
1922       /*output the fan of the grid points of the (i)th and (i-1)th grid line.
1923        */
1924       if(gridChain->getUlineIndex(i) < gridChain->getUlineIndex(i-1))
1925         {
1926           pStream->begin();
1927           pStream->insert(gridChain->get_vertex(i-1));
1928           for(j=gridChain->getUlineIndex(i); j<= gridChain->getUlineIndex(i-1); j++)
1929             pStream->insert(gridChain->getGrid()->get_u_value(j), gridChain->get_v_value(i));
1930           pStream->end(PRIMITIVE_STREAM_FAN);
1931         }
1932       else if(gridChain->getUlineIndex(i) > gridChain->getUlineIndex(i-1))
1933         {
1934           pStream->begin();
1935           pStream->insert(gridChain->get_vertex(i));
1936           for(j=gridChain->getUlineIndex(i); j>= gridChain->getUlineIndex(i-1); j--)
1937             pStream->insert(gridChain->getGrid()->get_u_value(j), gridChain->get_v_value(i-1));
1938           pStream->end(PRIMITIVE_STREAM_FAN);
1939         }
1940       /*otherwisem, the two are equal, so there is no fan to outout*/
1941     }
1942
1943   monoTriangulation2(upperVert, lowerVert, &vArray, 0, endIndex-beginIndex,
1944                      0, /*decreasing chain*/
1945                      pStream);
1946 }
1947
1948 /*return i, such that from begin to i-1 the chain is strictly u-monotone.
1949  */
1950 Int findIncreaseChainFromBegin(vertexArray* chain, Int begin ,Int end)
1951 {
1952   Int i=begin;
1953   Real prevU = chain->getVertex(i)[0];
1954   Real thisU;
1955   for(i=begin+1; i<=end; i++){
1956     thisU = chain->getVertex(i)[0];
1957
1958     if(prevU < thisU){
1959       prevU = thisU;
1960     }
1961     else
1962       break;
1963   }
1964   return i;
1965 }
1966
1967 /*check whether there is a vertex whose v value is strictly
1968  *inbetween vup vbelow
1969  *if no middle exists return -1, else return the idnex.
1970  */
1971 Int checkMiddle(vertexArray* chain, Int begin, Int end,
1972                 Real vup, Real vbelow)
1973 {
1974   Int i;
1975   for(i=begin; i<=end; i++)
1976     {
1977       if(chain->getVertex(i)[1] < vup && chain->getVertex(i)[1]>vbelow)
1978         return i;
1979     }
1980   return -1;
1981 }
1982
1983 /*the degenerat case of sampleLeftOneGridStep*/
1984 void sampleLeftOneGridStepNoMiddle(vertexArray* leftChain,
1985                                    Int beginLeftIndex,
1986                                    Int endLeftIndex,
1987                                    gridBoundaryChain* leftGridChain,
1988                                    Int leftGridChainStartIndex,
1989                                    primStream* pStream)
1990 {
1991   /*since there is no middle, there is at most one point which is on the
1992    *second grid line, there could be multiple points on the first (top)
1993    *grid line.
1994    */
1995
1996   leftGridChain->leftEndFan(leftGridChainStartIndex+1, pStream);
1997
1998   monoTriangulation2(leftGridChain->get_vertex(leftGridChainStartIndex),
1999                      leftGridChain->get_vertex(leftGridChainStartIndex+1),
2000                      leftChain,
2001                      beginLeftIndex,
2002                      endLeftIndex,
2003                      1, //is increase chain.
2004                      pStream);
2005 }
2006
2007
2008
2009 /*sampling the left area in between two grid lines.
2010  */
2011 void sampleLeftOneGridStep(vertexArray* leftChain,
2012                   Int beginLeftIndex,
2013                   Int endLeftIndex,
2014                   gridBoundaryChain* leftGridChain,
2015                   Int leftGridChainStartIndex,
2016                   primStream* pStream
2017                   )
2018 {
2019   if(checkMiddle(leftChain, beginLeftIndex, endLeftIndex,
2020                  leftGridChain->get_v_value(leftGridChainStartIndex),
2021                  leftGridChain->get_v_value(leftGridChainStartIndex+1))<0)
2022
2023     {
2024
2025       sampleLeftOneGridStepNoMiddle(leftChain, beginLeftIndex, endLeftIndex, leftGridChain, leftGridChainStartIndex, pStream);
2026       return;
2027     }
2028
2029   //copy into a polygon
2030   {
2031     directedLine* poly = NULL;
2032     sampledLine* sline;
2033     directedLine* dline;
2034     gridWrap* grid = leftGridChain->getGrid();
2035     Real vert1[2];
2036     Real vert2[2];
2037     Int i;
2038
2039     Int innerInd = leftGridChain->getInnerIndex(leftGridChainStartIndex+1);
2040     Int upperInd = leftGridChain->getUlineIndex(leftGridChainStartIndex);
2041     Int lowerInd = leftGridChain->getUlineIndex(leftGridChainStartIndex+1);
2042     Real upperV = leftGridChain->get_v_value(leftGridChainStartIndex);
2043     Real lowerV = leftGridChain->get_v_value(leftGridChainStartIndex+1);
2044
2045     //the upper gridline
2046     vert1[1] = vert2[1] = upperV;
2047     for(i=innerInd;     i>upperInd;     i--)
2048       {
2049         vert1[0]=grid->get_u_value(i);
2050         vert2[0]=grid->get_u_value(i-1);
2051         sline = new sampledLine(vert1, vert2);
2052         dline = new directedLine(INCREASING, sline);
2053         if(poly == NULL)
2054           poly = dline;
2055         else
2056           poly->insert(dline);
2057       }
2058
2059     //the edge connecting upper grid with left chain
2060     vert1[0] = grid->get_u_value(upperInd);
2061     vert1[1] = upperV;
2062     sline = new sampledLine(vert1, leftChain->getVertex(beginLeftIndex));
2063     dline = new directedLine(INCREASING, sline);
2064     if(poly == NULL)
2065       poly = dline;
2066     else
2067       poly->insert(dline);
2068
2069     //the left chain
2070     for(i=beginLeftIndex; i<endLeftIndex; i++)
2071       {
2072         sline = new sampledLine(leftChain->getVertex(i), leftChain->getVertex(i+1));
2073         dline = new directedLine(INCREASING, sline);
2074         poly->insert(dline);
2075       }
2076
2077     //the edge connecting left chain with lower gridline
2078     vert2[0] = grid->get_u_value(lowerInd);
2079     vert2[1] = lowerV;
2080     sline = new sampledLine(leftChain->getVertex(endLeftIndex), vert2);
2081     dline = new directedLine(INCREASING, sline);
2082     poly->insert(dline);
2083
2084     //the lower grid line
2085     vert1[1] = vert2[1] = lowerV;
2086     for(i=lowerInd; i<innerInd; i++)
2087       {
2088         vert1[0] = grid->get_u_value(i);
2089         vert2[0] = grid->get_u_value(i+1);
2090         sline = new sampledLine(vert1, vert2);
2091         dline = new directedLine(INCREASING, sline);
2092         poly->insert(dline);
2093       }
2094
2095     //the vertical grid line segement
2096     vert1[0]=vert2[0] = grid->get_u_value(innerInd);
2097     vert2[1]=upperV;
2098     vert1[1]=lowerV;
2099     sline=new sampledLine(vert1, vert2);
2100     dline=new directedLine(INCREASING, sline);
2101     poly->insert(dline);
2102     monoTriangulationOpt(poly, pStream);
2103     //cleanup
2104     poly->deleteSinglePolygonWithSline();
2105     return;
2106   }
2107
2108
2109
2110
2111
2112   Int i;
2113   if(1/*leftGridChain->getUlineIndex(leftGridChainStartIndex) >=
2114      leftGridChain->getUlineIndex(leftGridChainStartIndex+1)*/
2115      ) /*the second grid line is beyond the first one to the left*/
2116     {
2117       /*find the maximal U-monotone chain
2118        * of endLeftIndex, endLeftIndex-1, ...,
2119        */
2120       i=endLeftIndex;
2121       Real prevU = leftChain->getVertex(i)[0];
2122       for(i=endLeftIndex-1; i>=beginLeftIndex; i--){
2123         Real thisU = leftChain->getVertex(i)[0];
2124         if( prevU < thisU){
2125           prevU = thisU;
2126         }
2127         else
2128           break;
2129       }
2130       /*from endLeftIndex to i+1 is strictly U- monotone */
2131       /*if i+1==endLeftIndex and the vertex and leftchain is on the second gridline, then
2132        *we should use 2 vertices on the leftchain. If we only use one (endLeftIndex), then we
2133        *we would output degenerate triangles
2134        */
2135       if(i+1 == endLeftIndex && leftChain->getVertex(endLeftIndex)[1] == leftGridChain->get_v_value(1+leftGridChainStartIndex))
2136         i--;
2137
2138       Int j = beginLeftIndex/*endLeftIndex*/+1;
2139
2140
2141       if(leftGridChain->getInnerIndex(leftGridChainStartIndex+1) > leftGridChain->getUlineIndex(leftGridChainStartIndex))
2142         {
2143           j = findIncreaseChainFromBegin(leftChain, beginLeftIndex, i+1/*endLeftIndex*/);
2144
2145           Int temp = beginLeftIndex;
2146           /*now from begin to j-1 is strictly u-monotone*/
2147           /*if j-1 is on the first grid line, then we want to skip to the vertex which is strictly
2148            *below the grid line. This vertexmust exist since there is a 'corner turn' inbetween the two grid lines
2149            */
2150           if(j-1 == beginLeftIndex)
2151             {
2152               while(leftChain->getVertex(j-1)[1] == leftGridChain->get_v_value(leftGridChainStartIndex))
2153                 j++;
2154
2155               Real vert[2];
2156               vert[0] = leftGridChain->get_u_value(leftGridChainStartIndex);
2157               vert[1] = leftGridChain->get_v_value(leftGridChainStartIndex);
2158
2159               monoTriangulation2(
2160                                  vert/*leftChain->getVertex(beginLeftIndex)*/,
2161                                  leftChain->getVertex(j-1),
2162                                  leftChain,
2163                                  beginLeftIndex,
2164                                  j-2,
2165                                  1,
2166                                  pStream //increase chain
2167                                  );
2168
2169               temp = j-1;
2170             }
2171
2172           stripOfFanLeft(leftChain, j-1, temp/*beginLeftIndex*/, leftGridChain->getGrid(),
2173                          leftGridChain->getVlineIndex(leftGridChainStartIndex),
2174                          leftGridChain->getUlineIndex(leftGridChainStartIndex),
2175                          leftGridChain->getInnerIndex(leftGridChainStartIndex+1),
2176                          pStream,
2177                          1 /*the grid line is above the trim line*/
2178                          );
2179         }
2180
2181       stripOfFanLeft(leftChain, endLeftIndex, i+1, leftGridChain->getGrid(),
2182                      leftGridChain->getVlineIndex(leftGridChainStartIndex+1),
2183                      leftGridChain->getUlineIndex(leftGridChainStartIndex+1),
2184                      leftGridChain->getInnerIndex(leftGridChainStartIndex+1),
2185                      pStream,
2186                      0 /*the grid line is below the trim lines*/
2187                      );
2188
2189       /*monotone triangulate the remaining left chain togther with the
2190        *two vertices on the two grid v-lines.
2191        */
2192       Real vert[2][2];
2193       vert[0][0]=vert[1][0] = leftGridChain->getInner_u_value(leftGridChainStartIndex+1);
2194       vert[0][1] = leftGridChain->get_v_value(leftGridChainStartIndex);
2195       vert[1][1] = leftGridChain->get_v_value(leftGridChainStartIndex+1);
2196
2197 //      vertexArray right(vert, 2);
2198
2199       monoTriangulation2(
2200                          &vert[0][0], /*top vertex */
2201                          &vert[1][0], /*bottom vertex*/
2202                          leftChain,
2203                          /*beginLeftIndex*/j-1,
2204                          i+1,
2205                          1, /*an increasing chain*/
2206                          pStream);
2207     }
2208   else /*the second one is shorter than the first one to the left*/
2209     {
2210       /*find the maximal U-monotone chain of beginLeftIndex, beginLeftIndex+1,...,
2211        */
2212       i=beginLeftIndex;
2213       Real prevU = leftChain->getVertex(i)[0];
2214       for(i=beginLeftIndex+1; i<=endLeftIndex; i++){
2215         Real thisU = leftChain->getVertex(i)[0];
2216
2217         if(prevU < thisU){
2218           prevU = thisU;
2219         }
2220         else
2221           break;
2222       }
2223       /*from beginLeftIndex to i-1 is strictly U-monotone*/
2224
2225
2226       stripOfFanLeft(leftChain, i-1, beginLeftIndex, leftGridChain->getGrid(),
2227                          leftGridChain->getVlineIndex(leftGridChainStartIndex),
2228                          leftGridChain->getUlineIndex(leftGridChainStartIndex),
2229                          leftGridChain->getUlineIndex(leftGridChainStartIndex+1),
2230                          pStream,
2231                          1 /*the grid line is above the trim lines*/
2232                          );
2233       /*monotone triangulate the remaining left chain together with the
2234        *two vertices on the two grid v-lines.
2235        */
2236       Real vert[2][2];
2237       vert[0][0]=vert[1][0] = leftGridChain->get_u_value(leftGridChainStartIndex+1);
2238       vert[0][1] = leftGridChain->get_v_value(leftGridChainStartIndex);
2239       vert[1][1] = leftGridChain->get_v_value(leftGridChainStartIndex+1);
2240
2241       vertexArray right(vert, 2);
2242
2243       monoTriangulation2(
2244                          &vert[0][0], //top vertex
2245                          &vert[1][0], //bottom vertex
2246                          leftChain,
2247                          i-1,
2248                          endLeftIndex,
2249                          1, /*an increase chain*/
2250                          pStream);
2251
2252     }
2253 }
2254
2255 /*n_upper>=1
2256  *n_lower>=1
2257  */
2258 void triangulateXYMono(Int n_upper, Real upperVerts[][2],
2259                        Int n_lower, Real lowerVerts[][2],
2260                        primStream* pStream)
2261 {
2262   Int i,j,k,l;
2263   Real* leftMostV;
2264
2265   assert(n_upper>=1 && n_lower>=1);
2266   if(upperVerts[0][0] <= lowerVerts[0][0])
2267     {
2268       i=1;
2269       j=0;
2270       leftMostV = upperVerts[0];
2271     }
2272   else
2273     {
2274       i=0;
2275       j=1;
2276       leftMostV = lowerVerts[0];
2277     }
2278
2279   while(1)
2280     {
2281       if(i >= n_upper) /*case1: no more in upper*/
2282         {
2283
2284           if(j<n_lower-1) /*at least two vertices in lower*/
2285             {
2286               pStream->begin();
2287               pStream->insert(leftMostV);
2288               while(j<n_lower){
2289                 pStream->insert(lowerVerts[j]);
2290                 j++;
2291               }
2292               pStream->end(PRIMITIVE_STREAM_FAN);
2293             }
2294
2295           break;
2296         }
2297       else if(j>= n_lower) /*case2: no more in lower*/
2298         {
2299
2300           if(i<n_upper-1) /*at least two vertices in upper*/
2301             {
2302               pStream->begin();
2303               pStream->insert(leftMostV);
2304
2305               for(k=n_upper-1; k>=i; k--)
2306                 pStream->insert(upperVerts[k]);
2307
2308               pStream->end(PRIMITIVE_STREAM_FAN);
2309             }
2310
2311           break;
2312         }
2313       else /* case3: neither is empty, plus the leftMostV, there is at least one triangle to output*/
2314         {
2315
2316           if(upperVerts[i][0] <= lowerVerts[j][0])
2317             {
2318               pStream->begin();
2319               pStream->insert(lowerVerts[j]); /*the origin of this fan*/
2320
2321               /*find the last k>=i such that
2322                *upperverts[k][0] <= lowerverts[j][0]
2323                */
2324               k=i;
2325               while(k<n_upper)
2326                 {
2327                   if(upperVerts[k][0] > lowerVerts[j][0])
2328                     break;
2329                   k++;
2330                 }
2331               k--;
2332               for(l=k; l>=i; l--)/*the reverse is for two-face lighting*/
2333                 {
2334                   pStream->insert(upperVerts[l]);
2335                 }
2336               pStream->insert(leftMostV);
2337
2338               pStream->end(PRIMITIVE_STREAM_FAN);
2339               //update i for next loop
2340               i = k+1;
2341               leftMostV = upperVerts[k];
2342
2343             }
2344           else /*upperVerts[i][0] > lowerVerts[j][0]*/
2345             {
2346               pStream->begin();
2347               pStream->insert(upperVerts[i]);/*the origion of this fan*/
2348               pStream->insert(leftMostV);
2349               /*find the last k>=j such that
2350                *lowerverts[k][0] < upperverts[i][0]*/
2351               k=j;
2352               while(k< n_lower)
2353                 {
2354                   if(lowerVerts[k][0] >= upperVerts[i][0])
2355                     break;
2356                   pStream->insert(lowerVerts[k]);
2357                   k++;
2358                 }
2359               pStream->end(PRIMITIVE_STREAM_FAN);
2360               j=k;
2361               leftMostV = lowerVerts[j-1];
2362             }
2363         }
2364     }
2365 }
2366
2367
2368 void stripOfFanLeft(vertexArray* leftChain,
2369                     Int largeIndex,
2370                     Int smallIndex,
2371                     gridWrap* grid,
2372                     Int vlineIndex,
2373                     Int ulineSmallIndex,
2374                     Int ulineLargeIndex,
2375                     primStream* pStream,
2376                     Int gridLineUp /*1 if the grid line is above the trim lines*/
2377                      )
2378 {
2379   assert(largeIndex >= smallIndex);
2380
2381   Real grid_v_value;
2382   grid_v_value = grid->get_v_value(vlineIndex);
2383
2384   Real2* trimVerts=(Real2*) malloc(sizeof(Real2)* (largeIndex-smallIndex+1));
2385   assert(trimVerts);
2386
2387
2388   Real2* gridVerts=(Real2*) malloc(sizeof(Real2)* (ulineLargeIndex-ulineSmallIndex+1));
2389   assert(gridVerts);
2390
2391   Int k,i;
2392   if(gridLineUp) /*trim line is below grid line, so trim vertices are going right when index increases*/
2393     for(k=0, i=smallIndex; i<=largeIndex; i++, k++)
2394       {
2395       trimVerts[k][0] = leftChain->getVertex(i)[0];
2396       trimVerts[k][1] = leftChain->getVertex(i)[1];
2397     }
2398   else
2399     for(k=0, i=largeIndex; i>=smallIndex; i--, k++)
2400       {
2401         trimVerts[k][0] = leftChain->getVertex(i)[0];
2402         trimVerts[k][1] = leftChain->getVertex(i)[1];
2403       }
2404
2405   for(k=0, i=ulineSmallIndex; i<= ulineLargeIndex; i++, k++)
2406     {
2407       gridVerts[k][0] = grid->get_u_value(i);
2408       gridVerts[k][1] = grid_v_value;
2409     }
2410
2411   if(gridLineUp)
2412     triangulateXYMono(
2413                       ulineLargeIndex-ulineSmallIndex+1, gridVerts,
2414                       largeIndex-smallIndex+1, trimVerts,
2415                       pStream);
2416   else
2417     triangulateXYMono(largeIndex-smallIndex+1, trimVerts,
2418                       ulineLargeIndex-ulineSmallIndex+1, gridVerts,
2419                       pStream);
2420   free(trimVerts);
2421   free(gridVerts);
2422 }
2423
2424
2425
2426
2427