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