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: 2005/10/28 13:09:23 $ $Revision: 1.5 $
  35 */
  36 /*
  37 ** $Header: /home/krh/git/sync/mesa-cvs-repo/Mesa/src/glu/sgi/libnurbs/nurbtess/sampleMonoPoly.cc,v 1.5 2005/10/28 13:09:23 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       Int 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 = 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       Int tempI;
 521       //skip those points which are equal to v. (avoid degeneratcy)
 522       for(tempI = index1; tempI <= leftChainEndIndex; tempI++)
 523         if(leftChain->getVertex(tempI)[1] < v)
 524           break;
 525       if(tempI > leftChainEndIndex)
 526         ret_rightCornerWhere = 1;
 527       else
 528         {
 529           Real tempMax = leftChain->getVertex(tempI)[0];
 530           for(i=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 = 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 = 0;
 980   Int n_ulines = grid->get_n_ulines();
 981   Real uMin = grid->get_u_min();
 982   Real uMax = grid->get_u_max();
 983   /*
 984   Real vMin = grid->get_v_min();
 985   Real vMax = grid->get_v_max();
 986   */
 987   Real slop = 0.0, uinterc;
 988
 989 #ifdef SHORTEN_GRID_LINE
 990   //uintercBuf stores all the interction u value for each grid line
 991   //notice that lastGridIndex<= firstGridIndex
 992   Real *uintercBuf = (Real *) malloc (sizeof(Real) * (firstGridIndex-lastGridIndex+1));
 993   assert(uintercBuf);
 994 #endif
 995
 996   /*initialization to make vtail bigger than grid->...*/
 997   directedLine* dLine = topEdge;
 998   Real vtail = grid->get_v_value(firstGridIndex) + 1.0;
 999   Real tempMaxU = grid->get_u_min();
1000
1001
1002   /*for each grid line*/
1003   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1004     {
1005
1006       Real grid_v_value = grid->get_v_value(i);
1007
1008       /*check whether this grid line is below the current trim edge.*/
1009       if(vtail > grid_v_value)
1010         {
1011           /*since the grid line is below the trim edge, we
1012            *find the trim edge which will contain the trim line
1013            */
1014           while( (vtail=dLine->tail()[1]) > grid_v_value){
1015
1016             tempMaxU = max(tempMaxU, dLine->tail()[0]);
1017             dLine = dLine -> getNext();
1018           }
1019
1020           if( fabs(dLine->head()[1] - vtail) < ZERO)
1021             isHoriz = 1;
1022           else
1023             {
1024               isHoriz = 0;
1025               slop = (dLine->head()[0] - dLine->tail()[0]) / (dLine->head()[1]-vtail);
1026             }
1027         }
1028
1029       if(isHoriz)
1030         {
1031           uinterc = max(dLine->head()[0], dLine->tail()[0]);
1032         }
1033       else
1034         {
1035           uinterc = slop * (grid_v_value - vtail) + dLine->tail()[0];
1036         }
1037
1038       tempMaxU = max(tempMaxU, uinterc);
1039
1040       if(uinterc < uMin && uinterc >= uMin - ZERO)
1041         uinterc = uMin;
1042       if(uinterc > uMax && uinterc <= uMax + ZERO)
1043         uinterc = uMax;
1044
1045 #ifdef SHORTEN_GRID_LINE
1046       uintercBuf[k] = uinterc;
1047 #endif
1048
1049       assert(uinterc >= uMin && uinterc <= uMax);
1050        if(uinterc == uMax)
1051          ret_indices[k] = n_ulines-1;
1052        else
1053          ret_indices[k] = (Int)(((uinterc-uMin)/(uMax - uMin)) * (n_ulines-1)) + 1;
1054       if(ret_indices[k] >= n_ulines)
1055         ret_indices[k] = n_ulines-1;
1056
1057
1058       ret_innerIndices[k] = (Int)(((tempMaxU-uMin)/(uMax - uMin)) * (n_ulines-1)) + 1;
1059
1060       /*reinitialize tempMaxU for next grdiLine*/
1061       tempMaxU = uinterc;
1062     }
1063 #ifdef SHORTEN_GRID_LINE
1064   //for each grid line, compare the left grid point with the
1065   //intersection point. If the two points are too close, then
1066   //we should move the grid point one grid to the right
1067   //and accordingly we should update the inner index.
1068   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1069     {
1070       //check gridLine i
1071       //check ret_indices[k]
1072       Real a = grid->get_u_value(ret_indices[k]-1);
1073       Real b = grid->get_u_value(ret_indices[k]);
1074       assert(uintercBuf[k] >= a && uintercBuf < b);
1075       if( (b-uintercBuf[k]) <= 0.2 * (b-a)) //interc is very close to b
1076         {
1077           ret_indices[k]++;
1078         }
1079
1080       //check ret_innerIndices[k]
1081       if(k>0)
1082         {
1083           if(ret_innerIndices[k] < ret_indices[k-1])
1084             ret_innerIndices[k] = ret_indices[k-1];
1085           if(ret_innerIndices[k] < ret_indices[k])
1086             ret_innerIndices[k] = ret_indices[k];
1087         }
1088     }
1089   //clean up
1090   free(uintercBuf);
1091 #endif
1092 }
1093
1094 void findRightGridIndices(directedLine* topEdge, Int firstGridIndex, Int lastGridIndex, gridWrap* grid,  Int* ret_indices, Int* ret_innerIndices)
1095 {
1096
1097   Int i,k;
1098   Int n_ulines = grid->get_n_ulines();
1099   Real uMin = grid->get_u_min();
1100   Real uMax = grid->get_u_max();
1101   /*
1102   Real vMin = grid->get_v_min();
1103   Real vMax = grid->get_v_max();
1104   */
1105   Real slop = 0.0, uinterc;
1106
1107 #ifdef SHORTEN_GRID_LINE
1108   //uintercBuf stores all the interction u value for each grid line
1109   //notice that firstGridIndex >= lastGridIndex
1110   Real *uintercBuf = (Real *) malloc (sizeof(Real) * (firstGridIndex-lastGridIndex+1));
1111   assert(uintercBuf);
1112 #endif
1113
1114   /*initialization to make vhead bigger than grid->v_value...*/
1115   directedLine* dLine = topEdge->getPrev();
1116   Real vhead = dLine->tail()[1];
1117   Real tempMinU = grid->get_u_max();
1118
1119   /*for each grid line*/
1120   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1121     {
1122
1123       Real grid_v_value = grid->get_v_value(i);
1124
1125
1126       /*check whether this grid line is below the current trim edge.*/
1127       if(vhead >= grid_v_value)
1128         {
1129           /*since the grid line is below the tail of the trim edge, we
1130            *find the trim edge which will contain the trim line
1131            */
1132           while( (vhead=dLine->head()[1]) > grid_v_value){
1133             tempMinU = min(tempMinU, dLine->head()[0]);
1134             dLine = dLine -> getPrev();
1135           }
1136
1137           /*skip the equality in the case of degenerat case: horizontal */
1138           while(dLine->head()[1] == grid_v_value)
1139             dLine = dLine->getPrev();
1140
1141           assert( dLine->tail()[1] != dLine->head()[1]);
1142           slop = (dLine->tail()[0] - dLine->head()[0]) / (dLine->tail()[1]-dLine->head()[1]);
1143           /*
1144            if(dLine->tail()[1] == vhead)
1145              isHoriz = 1;
1146              else
1147             {
1148               isHoriz = 0;
1149               slop = (dLine->tail()[0] - dLine->head()[0]) / (dLine->tail()[1]-vhead);
1150             }
1151             */
1152         }
1153         uinterc = slop * (grid_v_value - dLine->head()[1]) + dLine->head()[0];
1154
1155       //in case unterc is outside of the grid due to floating point
1156       if(uinterc < uMin)
1157         uinterc = uMin;
1158       else if(uinterc > uMax)
1159         uinterc = uMax;
1160
1161 #ifdef SHORTEN_GRID_LINE
1162       uintercBuf[k] = uinterc;
1163 #endif
1164
1165       tempMinU = min(tempMinU, uinterc);
1166
1167       assert(uinterc >= uMin && uinterc <= uMax);
1168
1169       if(uinterc == uMin)
1170         ret_indices[k] = 0;
1171       else
1172         ret_indices[k] = (int)ceil((((uinterc-uMin)/(uMax - uMin)) * (n_ulines-1))) -1;
1173 /*
1174 if(ret_indices[k] >= grid->get_n_ulines())
1175   {
1176   printf("ERROR3\n");
1177   exit(0);
1178 }
1179 if(ret_indices[k] < 0)
1180   {
1181   printf("ERROR4\n");
1182   exit(0);
1183 }
1184 */
1185       ret_innerIndices[k] = (int)ceil ((((tempMinU-uMin)/(uMax - uMin)) * (n_ulines-1))) -1;
1186
1187       tempMinU = uinterc;
1188     }
1189 #ifdef SHORTEN_GRID_LINE
1190   //for each grid line, compare the left grid point with the
1191   //intersection point. If the two points are too close, then
1192   //we should move the grid point one grid to the right
1193   //and accordingly we should update the inner index.
1194   for(k=0, i=firstGridIndex; i>=lastGridIndex; i--, k++)
1195     {
1196       //check gridLine i
1197       //check ret_indices[k]
1198       Real a = grid->get_u_value(ret_indices[k]);
1199       Real b = grid->get_u_value(ret_indices[k]+1);
1200       assert(uintercBuf[k] > a && uintercBuf <= b);
1201       if( (uintercBuf[k]-a) <= 0.2 * (b-a)) //interc is very close to a
1202         {
1203           ret_indices[k]--;
1204         }
1205
1206       //check ret_innerIndices[k]
1207       if(k>0)
1208         {
1209           if(ret_innerIndices[k] > ret_indices[k-1])
1210             ret_innerIndices[k] = ret_indices[k-1];
1211           if(ret_innerIndices[k] > ret_indices[k])
1212             ret_innerIndices[k] = ret_indices[k];
1213         }
1214     }
1215   //clean up
1216   free(uintercBuf);
1217 #endif
1218 }
1219
1220
1221 void sampleMonoPoly(directedLine* polygon, gridWrap* grid, Int ulinear, Int vlinear, primStream* pStream, rectBlockArray* rbArray)
1222 {
1223 /*
1224 {
1225 grid->print();
1226 polygon->writeAllPolygons("zloutputFile");
1227 exit(0);
1228 }
1229 */
1230
1231 if(grid->get_n_ulines() == 2 ||
1232    grid->get_n_vlines() == 2)
1233 {
1234   if(ulinear && grid->get_n_ulines() == 2)
1235     {
1236       monoTriangulationFun(polygon, compV2InY, pStream);
1237       return;
1238     }
1239   else if(DBG_isConvex(polygon) && polygon->numEdges() >=4)
1240      {
1241        triangulateConvexPoly(polygon, ulinear, vlinear, pStream);
1242        return;
1243      }
1244   else if(vlinear || DBG_is_U_direction(polygon))
1245     {
1246       Int n_cusps;//num interior cusps
1247       Int n_edges = polygon->numEdges();
1248       directedLine** cusps = (directedLine**) malloc(sizeof(directedLine*) * n_edges);
1249       assert(cusps);
1250       findInteriorCuspsX(polygon, n_cusps, cusps);
1251
1252       if(n_cusps == 0) //u_monotone
1253         {
1254
1255           monoTriangulationFun(polygon, compV2InX, pStream);
1256
1257           free(cusps);
1258           return;
1259         }
1260       else if(n_cusps == 1) //one interior cusp
1261         {
1262
1263           directedLine* new_polygon = polygonConvert(cusps[0]);
1264
1265           directedLine* other = findDiagonal_singleCuspX( new_polygon);
1266
1267
1268
1269           //<other> should NOT be null unless there are self-intersecting
1270           //trim curves. In that case, we don't want to core dump, instead,
1271           //we triangulate anyway, and print out error message.
1272           if(other == NULL)
1273             {
1274               monoTriangulationFun(polygon, compV2InX, pStream);
1275               free(cusps);
1276               return;
1277             }
1278
1279           directedLine* ret_p1;
1280           directedLine* ret_p2;
1281
1282           new_polygon->connectDiagonal_2slines(new_polygon, other,
1283                                            &ret_p1,
1284                                            &ret_p2,
1285                                            new_polygon);
1286
1287           monoTriangulationFun(ret_p1, compV2InX, pStream);
1288           monoTriangulationFun(ret_p2, compV2InX, pStream);
1289
1290           ret_p1->deleteSinglePolygonWithSline();
1291           ret_p2->deleteSinglePolygonWithSline();
1292
1293           free(cusps);
1294           return;
1295          }
1296      free(cusps);
1297      }
1298 }
1299
1300   /*find the top and bottom of the polygon. It is supposed to be
1301    *a V-monotone polygon
1302    */
1303
1304   directedLine* tempV;
1305   directedLine* topV;
1306   directedLine* botV;
1307   topV = botV = polygon;
1308
1309   for(tempV = polygon->getNext(); tempV != polygon; tempV = tempV->getNext())
1310     {
1311       if(compV2InY(topV->head(), tempV->head())<0) {
1312
1313         topV = tempV;
1314       }
1315       if(compV2InY(botV->head(), tempV->head())>0) {
1316
1317         botV = tempV;
1318       }
1319     }
1320
1321   /*find the first(top) and the last (bottom) grid line which intersect the
1322    *this polygon
1323    */
1324   Int firstGridIndex; /*the index in the grid*/
1325   Int lastGridIndex;
1326   firstGridIndex = (Int) ((topV->head()[1] - grid->get_v_min()) / (grid->get_v_max() - grid->get_v_min()) * (grid->get_n_vlines()-1));
1327   lastGridIndex  = (Int) ((botV->head()[1] - grid->get_v_min()) / (grid->get_v_max() - grid->get_v_min()) * (grid->get_n_vlines()-1)) + 1;
1328
1329
1330   /*find the interval inside  the polygon for each gridline*/
1331   Int *leftGridIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1332   assert(leftGridIndices);
1333   Int *rightGridIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1334   assert(rightGridIndices);
1335   Int *leftGridInnerIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1336   assert(leftGridInnerIndices);
1337   Int *rightGridInnerIndices = (Int*) malloc(sizeof(Int) * (firstGridIndex - lastGridIndex +1));
1338   assert(rightGridInnerIndices);
1339
1340   findLeftGridIndices(topV, firstGridIndex, lastGridIndex, grid,  leftGridIndices, leftGridInnerIndices);
1341
1342   findRightGridIndices(topV, firstGridIndex, lastGridIndex, grid,  rightGridIndices, rightGridInnerIndices);
1343
1344   gridBoundaryChain leftGridChain(grid, firstGridIndex, firstGridIndex-lastGridIndex+1, leftGridIndices, leftGridInnerIndices);
1345
1346   gridBoundaryChain rightGridChain(grid, firstGridIndex, firstGridIndex-lastGridIndex+1, rightGridIndices, rightGridInnerIndices);
1347
1348
1349
1350 //  leftGridChain.draw();
1351 //  leftGridChain.drawInner();
1352 //  rightGridChain.draw();
1353 //  rightGridChain.drawInner();
1354   /*(1) determine the grid boundaries (left and right).
1355    *(2) process polygon  into two monotone chaines: use vertexArray
1356    *(3) call sampleMonoPolyRec
1357    */
1358
1359   /*copy the two chains into vertexArray datastructure*/
1360   Int i;
1361   vertexArray leftChain(20); /*this is a dynamic array*/
1362   for(i=1; i<=topV->get_npoints()-2; i++) { /*the first vertex is the top vertex which doesn't belong to inc_chain*/
1363     leftChain.appendVertex(topV->getVertex(i));
1364   }
1365   for(tempV = topV->getNext(); tempV != botV; tempV = tempV->getNext())
1366     {
1367       for(i=0; i<=tempV->get_npoints()-2; i++){
1368         leftChain.appendVertex(tempV->getVertex(i));
1369       }
1370     }
1371
1372   vertexArray rightChain(20);
1373   for(tempV = topV->getPrev(); tempV != botV; tempV = tempV->getPrev())
1374     {
1375       for(i=tempV->get_npoints()-2; i>=0; i--){
1376         rightChain.appendVertex(tempV->getVertex(i));
1377       }
1378     }
1379   for(i=botV->get_npoints()-2; i>=1; i--){
1380     rightChain.appendVertex(tempV->getVertex(i));
1381   }
1382
1383   sampleMonoPolyRec(topV->head(),
1384                     botV->head(),
1385                     &leftChain,
1386                     0,
1387                     &rightChain,
1388                     0,
1389                     &leftGridChain,
1390                     &rightGridChain,
1391                     0,
1392                     pStream,
1393                     rbArray);
1394
1395
1396   /*cleanup space*/
1397   free(leftGridIndices);
1398   free(rightGridIndices);
1399   free(leftGridInnerIndices);
1400   free(rightGridInnerIndices);
1401 }
1402
1403 void sampleMonoPolyRec(
1404                        Real* topVertex,
1405                        Real* botVertex,
1406                        vertexArray* leftChain,
1407                        Int leftStartIndex,
1408                        vertexArray* rightChain,
1409                        Int rightStartIndex,
1410                        gridBoundaryChain* leftGridChain,
1411                        gridBoundaryChain* rightGridChain,
1412                        Int gridStartIndex,
1413                        primStream* pStream,
1414                        rectBlockArray* rbArray)
1415 {
1416
1417   /*find the first connected component, and the four corners.
1418    */
1419   Int index1, index2; /*the first and last grid line of the first connected component*/
1420
1421   if(topVertex[1] <= botVertex[1])
1422     return;
1423
1424   /*find i so that the grid line is below the top vertex*/
1425   Int i=gridStartIndex;
1426   while (i < leftGridChain->get_nVlines())
1427     {
1428       if(leftGridChain->get_v_value(i) < topVertex[1])
1429         break;
1430       i++;
1431     }
1432
1433   /*find the first connected component starting with i*/
1434   /*find index1 so that left_uline_index <= right_uline_index, that is, this
1435    *grid line contains at least one inner grid point
1436    */
1437   index1=i;
1438   int num_skipped_grid_lines=0;
1439   while(index1 < leftGridChain->get_nVlines())
1440     {
1441       if(leftGridChain->getUlineIndex(index1) <= rightGridChain->getUlineIndex(index1))
1442         break;
1443       num_skipped_grid_lines++;
1444       index1++;
1445     }
1446
1447
1448
1449   if(index1 >= leftGridChain->get_nVlines()) /*no grid line exists which has inner point*/
1450     {
1451       /*stop recursion, ...*/
1452       /*monotone triangulate it...*/
1453 //      printf("no grid line exists\n");
1454 /*
1455       monoTriangulationRecOpt(topVertex, botVertex, leftChain, leftStartIndex,
1456                            rightChain, rightStartIndex, pStream);
1457 */
1458
1459 if(num_skipped_grid_lines <2)
1460   {
1461     monoTriangulationRecGenOpt(topVertex, botVertex, leftChain, leftStartIndex,
1462                                leftChain->getNumElements()-1,
1463                                rightChain, rightStartIndex,
1464                                rightChain->getNumElements()-1,
1465                                pStream);
1466   }
1467 else
1468   {
1469     //the optimum way to triangulate is top-down since this polygon
1470     //is narrow-long.
1471     monoTriangulationRec(topVertex, botVertex, leftChain, leftStartIndex,
1472                          rightChain, rightStartIndex, pStream);
1473   }
1474
1475 /*
1476       monoTriangulationRec(topVertex, botVertex, leftChain, leftStartIndex,
1477                            rightChain, rightStartIndex, pStream);
1478 */
1479
1480 /*      monoTriangulationRecGenTBOpt(topVertex, botVertex,
1481                                    leftChain, leftStartIndex, leftChain->getNumElements()-1,
1482                                    rightChain, rightStartIndex, rightChain->getNumElements()-1,
1483                                    pStream);*/
1484
1485
1486
1487     }
1488   else
1489     {
1490
1491       /*find index2 so that left_inner_index <= right_inner_index holds until index2*/
1492       index2=index1+1;
1493       if(index2 < leftGridChain->get_nVlines())
1494         while(leftGridChain->getInnerIndex(index2) <= rightGridChain->getInnerIndex(index2))
1495           {
1496             index2++;
1497             if(index2 >= leftGridChain->get_nVlines())
1498               break;
1499           }
1500
1501       index2--;
1502
1503
1504
1505       /*the neck*/
1506       Int neckLeftIndex;
1507       Int neckRightIndex;
1508
1509       /*the four corners*/
1510       Int up_leftCornerWhere;
1511       Int up_leftCornerIndex;
1512       Int up_rightCornerWhere;
1513       Int up_rightCornerIndex;
1514       Int down_leftCornerWhere;
1515       Int down_leftCornerIndex;
1516       Int down_rightCornerWhere;
1517       Int down_rightCornerIndex;
1518
1519       Real* tempBotVertex; /*the bottom vertex for this component*/
1520       Real* nextTopVertex=NULL; /*for the recursion*/
1521       Int nextLeftStartIndex=0;
1522       Int nextRightStartIndex=0;
1523
1524       /*find the points below the grid line index2 on both chains*/
1525       Int botLeftIndex = leftChain->findIndexStrictBelowGen(
1526                                                       leftGridChain->get_v_value(index2),
1527                                                       leftStartIndex,
1528                                                       leftChain->getNumElements()-1);
1529       Int botRightIndex = rightChain->findIndexStrictBelowGen(
1530                                                         rightGridChain->get_v_value(index2),
1531                                                         rightStartIndex,
1532                                                         rightChain->getNumElements()-1);
1533       /*if either botLeftIndex>= numelements,
1534        *        or botRightIndex >= numelemnet,
1535        *then there is no neck exists. the bottom vertex is botVertex,
1536        */
1537       if(! findNeckF(leftChain, botLeftIndex, rightChain, botRightIndex,
1538                    leftGridChain, rightGridChain, index2, neckLeftIndex, neckRightIndex))
1539          /*
1540       if(botLeftIndex == leftChain->getNumElements() ||
1541          botRightIndex == rightChain->getNumElements())
1542          */
1543         {
1544 #ifdef MYDEBUG
1545           printf("neck NOT exists, botRightIndex=%i\n", botRightIndex);
1546 #endif
1547
1548           tempBotVertex =  botVertex;
1549           nextTopVertex = botVertex;
1550           botLeftIndex = leftChain->getNumElements()-1;
1551           botRightIndex = rightChain->getNumElements()-1;
1552         }
1553       else /*neck exists*/
1554         {
1555 #ifdef MYDEBUG
1556           printf("neck exists\n");
1557 #endif
1558
1559           /*
1560           findNeck(leftChain, botLeftIndex,
1561                    rightChain, botRightIndex,
1562                    neckLeftIndex,
1563                    neckRightIndex);
1564                    */
1565 #ifdef MYDEBUG
1566 printf("neck is found, neckLeftIndex=%i, neckRightIndex=%i\n", neckLeftIndex, neckRightIndex);
1567 glBegin(GL_LINES);
1568 glVertex2fv(leftChain->getVertex(neckLeftIndex));
1569 glVertex2fv(rightChain->getVertex(neckRightIndex));
1570 glEnd();
1571 #endif
1572
1573           if(leftChain->getVertex(neckLeftIndex)[1] <= rightChain->getVertex(neckRightIndex)[1])
1574             {
1575               tempBotVertex = leftChain->getVertex(neckLeftIndex);
1576               botLeftIndex = neckLeftIndex-1;
1577               botRightIndex = neckRightIndex;
1578               nextTopVertex = rightChain->getVertex(neckRightIndex);
1579               nextLeftStartIndex = neckLeftIndex;
1580               nextRightStartIndex = neckRightIndex+1;
1581             }
1582           else
1583             {
1584               tempBotVertex = rightChain->getVertex(neckRightIndex);
1585               botLeftIndex = neckLeftIndex;
1586               botRightIndex = neckRightIndex-1;
1587               nextTopVertex = leftChain->getVertex(neckLeftIndex);
1588               nextLeftStartIndex = neckLeftIndex+1;
1589               nextRightStartIndex = neckRightIndex;
1590             }
1591         }
1592
1593       findUpCorners(topVertex,
1594                     leftChain,
1595                     leftStartIndex, botLeftIndex,
1596                     rightChain,
1597                     rightStartIndex, botRightIndex,
1598                     leftGridChain->get_v_value(index1),
1599                     leftGridChain->get_u_value(index1),
1600                     rightGridChain->get_u_value(index1),
1601                     up_leftCornerWhere,
1602                     up_leftCornerIndex,
1603                     up_rightCornerWhere,
1604                     up_rightCornerIndex);
1605
1606       findDownCorners(tempBotVertex,
1607                       leftChain,
1608                       leftStartIndex, botLeftIndex,
1609                       rightChain,
1610                       rightStartIndex, botRightIndex,
1611                       leftGridChain->get_v_value(index2),
1612                       leftGridChain->get_u_value(index2),
1613                       rightGridChain->get_u_value(index2),
1614                       down_leftCornerWhere,
1615                       down_leftCornerIndex,
1616                       down_rightCornerWhere,
1617                       down_rightCornerIndex);
1618 #ifdef MYDEBUG
1619       printf("find corners done, down_leftwhere=%i, down_righwhere=%i,\n",down_leftCornerWhere, down_rightCornerWhere );
1620       printf("find corners done, up_leftwhere=%i, up_righwhere=%i,\n",up_leftCornerWhere, up_rightCornerWhere );
1621       printf("find corners done, up_leftindex=%i, up_righindex=%i,\n",up_leftCornerIndex, up_rightCornerIndex );
1622       printf("find corners done, down_leftindex=%i, down_righindex=%i,\n",down_leftCornerIndex, down_rightCornerIndex );
1623 #endif
1624
1625 /*
1626       drawCorners(topVertex,
1627                   tempBotVertex,
1628                   leftChain,
1629                   rightChain,
1630                   leftGridChain,
1631                   rightGridChain,
1632                   index1,
1633                   index2,
1634                   up_leftCornerWhere,
1635                   up_leftCornerIndex,
1636                   up_rightCornerWhere,
1637                   up_rightCornerIndex,
1638                   down_leftCornerWhere,
1639                   down_leftCornerIndex,
1640                   down_rightCornerWhere,
1641                   down_rightCornerIndex);
1642 */
1643
1644
1645       sampleConnectedComp(topVertex, tempBotVertex,
1646                           leftChain,
1647                           leftStartIndex, botLeftIndex,
1648                           rightChain,
1649                           rightStartIndex, botRightIndex,
1650                           leftGridChain,
1651                           rightGridChain,
1652                           index1, index2,
1653                           up_leftCornerWhere,
1654                           up_leftCornerIndex,
1655                           up_rightCornerWhere,
1656                           up_rightCornerIndex,
1657                           down_leftCornerWhere,
1658                           down_leftCornerIndex,
1659                           down_rightCornerWhere,
1660                           down_rightCornerIndex,
1661                           pStream,
1662                           rbArray
1663                           );
1664
1665       /*recursion*/
1666
1667       sampleMonoPolyRec(
1668                         nextTopVertex,
1669                         botVertex,
1670                         leftChain,
1671                         nextLeftStartIndex,
1672                         rightChain,
1673                         nextRightStartIndex,
1674                         leftGridChain,
1675                         rightGridChain,
1676                         index2+1,
1677                         pStream, rbArray);
1678
1679
1680     }
1681
1682 }
1683
1684 void sampleLeftStrip(vertexArray* leftChain,
1685                      Int topLeftIndex,
1686                      Int botLeftIndex,
1687                      gridBoundaryChain* leftGridChain,
1688                      Int leftGridChainStartIndex,
1689                      Int leftGridChainEndIndex,
1690                      primStream* pStream
1691                      )
1692 {
1693   assert(leftChain->getVertex(topLeftIndex)[1] > leftGridChain->get_v_value(leftGridChainStartIndex));
1694   assert(leftChain->getVertex(topLeftIndex+1)[1] <= leftGridChain->get_v_value(leftGridChainStartIndex));
1695   assert(leftChain->getVertex(botLeftIndex)[1] <= leftGridChain->get_v_value(leftGridChainEndIndex));
1696   assert(leftChain->getVertex(botLeftIndex-1)[1] > leftGridChain->get_v_value(leftGridChainEndIndex));
1697
1698   /*
1699    *(1)find the last grid line which doesn'; pass below
1700    * this first edge, sample this region: one trim edge and
1701    * possily multiple grid lines.
1702    */
1703   Real *upperVert, *lowerVert; /*the end points of the first trim edge*/
1704   upperVert = leftChain->getVertex(topLeftIndex);
1705   lowerVert = leftChain->getVertex(topLeftIndex+1);
1706
1707   Int index = leftGridChainStartIndex;
1708   while(leftGridChain->get_v_value(index) >= lowerVert[1]){
1709     index++;
1710     if(index > leftGridChainEndIndex)
1711       break;
1712   }
1713   index--;
1714
1715   sampleLeftSingleTrimEdgeRegion(upperVert, lowerVert,
1716                                  leftGridChain,
1717                                  leftGridChainStartIndex,
1718                                  index,
1719                                  pStream);
1720   sampleLeftStripRec(leftChain, topLeftIndex+1, botLeftIndex,
1721                      leftGridChain, index, leftGridChainEndIndex,
1722                      pStream);
1723
1724 }
1725
1726 void sampleLeftStripRec(vertexArray* leftChain,
1727                      Int topLeftIndex,
1728                      Int botLeftIndex,
1729                      gridBoundaryChain* leftGridChain,
1730                      Int leftGridChainStartIndex,
1731                      Int leftGridChainEndIndex,
1732                         primStream* pStream
1733                      )
1734 {
1735   /*now top left trim vertex is below the top grid line.
1736    */
1737   /*stop condition: if topLeftIndex >= botLeftIndex, then stop.
1738    */
1739   if(topLeftIndex >= botLeftIndex)
1740     return;
1741
1742   /*find the last trim vertex which is above the second top grid line:
1743    * index1.
1744    *and sampleLeftOneGridStep(leftchain, topLeftIndex, index1, leftGridChain,
1745    * leftGridChainStartIndex).
1746    * index1 could be equal to topLeftIndex.
1747    */
1748   Real secondGridChainV = leftGridChain->get_v_value(leftGridChainStartIndex+1);
1749   assert(leftGridChainStartIndex < leftGridChainEndIndex);
1750   Int index1 = topLeftIndex;
1751   while(leftChain->getVertex(index1)[1] > secondGridChainV)
1752     index1++;
1753   index1--;
1754
1755   sampleLeftOneGridStep(leftChain, topLeftIndex, index1, leftGridChain, leftGridChainStartIndex, pStream);
1756
1757
1758   /*
1759    * Let the next trim vertex be nextTrimVertIndex (which should be
1760    *  below the second grid line).
1761    * Find the last grid line index2 which is above nextTrimVert.
1762    * sampleLeftSingleTrimEdgeRegion(uppervert[2], lowervert[2],
1763    *                      leftGridChain, leftGridChainStartIndex+1, index2).
1764    */
1765   Real *uppervert, *lowervert;
1766   uppervert = leftChain->getVertex(index1);
1767   lowervert = leftChain->getVertex(index1+1);
1768   Int index2 = leftGridChainStartIndex+1;
1769
1770   while(leftGridChain->get_v_value(index2) >= lowervert[1])
1771     {
1772       index2++;
1773       if(index2 > leftGridChainEndIndex)
1774         break;
1775     }
1776   index2--;
1777   sampleLeftSingleTrimEdgeRegion(uppervert, lowervert, leftGridChain, leftGridChainStartIndex+1, index2,  pStream);
1778
1779    /* sampleLeftStripRec(leftChain,
1780                         nextTrimVertIndex,
1781                         botLeftIndex,
1782                         leftGridChain,
1783                         index2,
1784                         leftGridChainEndIndex
1785                         )
1786    *
1787    */
1788   sampleLeftStripRec(leftChain, index1+1, botLeftIndex, leftGridChain, index2, leftGridChainEndIndex, pStream);
1789
1790 }
1791
1792
1793 /***************begin RecF***********************/
1794 /* the gridlines from leftGridChainStartIndex to
1795  * leftGridChainEndIndex are assumed to form a
1796  * connected component.
1797  * the trim vertex of topLeftIndex is assumed to
1798  * be below the first gridline, and the tim vertex
1799  * of botLeftIndex is assumed to be above the last
1800  * grid line.
1801  * If botLeftIndex < topLeftIndex, then no connected componeent exists, and this funcion returns without
1802  * outputing any triangles.
1803  * Otherwise botLeftIndex >= topLeftIndex, there is at least one triangle to output.
1804  */
1805 void sampleLeftStripRecF(vertexArray* leftChain,
1806                      Int topLeftIndex,
1807                      Int botLeftIndex,
1808                      gridBoundaryChain* leftGridChain,
1809                      Int leftGridChainStartIndex,
1810                      Int leftGridChainEndIndex,
1811                         primStream* pStream
1812                      )
1813 {
1814   /*now top left trim vertex is below the top grid line.
1815    */
1816   /*stop condition: if topLeftIndex > botLeftIndex, then stop.
1817    */
1818   if(topLeftIndex > botLeftIndex)
1819     return;
1820
1821   /*if there is only one grid Line, return.*/
1822
1823   if(leftGridChainStartIndex>=leftGridChainEndIndex)
1824     return;
1825
1826
1827   assert(leftChain->getVertex(topLeftIndex)[1] <= leftGridChain->get_v_value(leftGridChainStartIndex) &&
1828          leftChain->getVertex(botLeftIndex)[1] >= leftGridChain->get_v_value(leftGridChainEndIndex));
1829
1830   /*firs find the first trim vertex which is below or equal to the second top grid line:
1831    * index1.
1832    */
1833   Real secondGridChainV = leftGridChain->get_v_value(leftGridChainStartIndex+1);
1834
1835
1836   Int index1 = topLeftIndex;
1837
1838   while(leftChain->getVertex(index1)[1] > secondGridChainV){
1839     index1++;
1840     if(index1>botLeftIndex)
1841       break;
1842   }
1843
1844   /*now leftChain->getVertex(index-1)[1] > secondGridChainV and
1845    *    leftChain->getVertex(index)[1] <= secondGridChainV
1846    *If equality holds, then we should include the vertex index1, otherwise we include only index1-1, to
1847    *perform sampleOneGridStep.
1848    */
1849   if(index1>botLeftIndex)
1850     index1--;
1851   else if(leftChain->getVertex(index1)[1] < secondGridChainV)
1852     index1--;
1853
1854   /*now we have leftChain->getVertex(index1)[1] >= secondGridChainV, and
1855    *  leftChain->getVertex(index1+1)[1] <= secondGridChainV
1856    */
1857
1858
1859   sampleLeftOneGridStep(leftChain, topLeftIndex, index1, leftGridChain, leftGridChainStartIndex, pStream);
1860
1861
1862   /*if leftChain->getVertex(index1)[1] == secondGridChainV, then we can recursively do the rest.
1863    */
1864   if(leftChain->getVertex(index1)[1] == secondGridChainV)
1865     {
1866
1867       sampleLeftStripRecF(leftChain, index1, botLeftIndex,leftGridChain, leftGridChainStartIndex+1, leftGridChainEndIndex, pStream);
1868     }
1869   else if(index1 < botLeftIndex)
1870     {
1871
1872       /* Otherwise, we have leftChain->getVertex(index1)[1] > secondGridChainV,
1873        * let the next trim vertex be nextTrimVertIndex (which should be  strictly
1874        *  below the second grid line).
1875        * Find the last grid line index2 which is above nextTrimVert.
1876        * sampleLeftSingleTrimEdgeRegion(uppervert[2], lowervert[2],
1877        *                      leftGridChain, leftGridChainStartIndex+1, index2).
1878        */
1879       Real *uppervert, *lowervert;
1880       uppervert = leftChain->getVertex(index1);
1881       lowervert = leftChain->getVertex(index1+1); //okay since index1<botLeftIndex
1882       Int index2 = leftGridChainStartIndex+1;
1883
1884
1885       while(leftGridChain->get_v_value(index2) >= lowervert[1])
1886         {
1887           index2++;
1888           if(index2 > leftGridChainEndIndex)
1889             break;
1890         }
1891       index2--;
1892
1893
1894       sampleLeftSingleTrimEdgeRegion(uppervert, lowervert, leftGridChain, leftGridChainStartIndex+1, index2,  pStream);
1895
1896       /*recursion*/
1897
1898       sampleLeftStripRecF(leftChain, index1+1, botLeftIndex, leftGridChain, index2, leftGridChainEndIndex, pStream);
1899     }
1900
1901 }
1902
1903 /***************End RecF***********************/
1904
1905 /*sample the left area in between one trim edge and multiple grid lines.
1906  * all the grid lines should be in between the two end poins of the
1907  *trim edge.
1908  */
1909 void sampleLeftSingleTrimEdgeRegion(Real upperVert[2], Real lowerVert[2],
1910                                     gridBoundaryChain* gridChain,
1911                                     Int beginIndex,
1912                                     Int endIndex,
1913                                     primStream* pStream)
1914 {
1915   Int i,j,k;
1916
1917   vertexArray vArray(endIndex-beginIndex+1);
1918   vArray.appendVertex(gridChain->get_vertex(beginIndex));
1919
1920   for(k=1, i=beginIndex+1; i<=endIndex; i++, k++)
1921     {
1922       vArray.appendVertex(gridChain->get_vertex(i));
1923
1924       /*output the fan of the grid points of the (i)th and (i-1)th grid line.
1925        */
1926       if(gridChain->getUlineIndex(i) < gridChain->getUlineIndex(i-1))
1927         {
1928           pStream->begin();
1929           pStream->insert(gridChain->get_vertex(i-1));
1930           for(j=gridChain->getUlineIndex(i); j<= gridChain->getUlineIndex(i-1); j++)
1931             pStream->insert(gridChain->getGrid()->get_u_value(j), gridChain->get_v_value(i));
1932           pStream->end(PRIMITIVE_STREAM_FAN);
1933         }
1934       else if(gridChain->getUlineIndex(i) > gridChain->getUlineIndex(i-1))
1935         {
1936           pStream->begin();
1937           pStream->insert(gridChain->get_vertex(i));
1938           for(j=gridChain->getUlineIndex(i); j>= gridChain->getUlineIndex(i-1); j--)
1939             pStream->insert(gridChain->getGrid()->get_u_value(j), gridChain->get_v_value(i-1));
1940           pStream->end(PRIMITIVE_STREAM_FAN);
1941         }
1942       /*otherwisem, the two are equal, so there is no fan to outout*/
1943     }
1944
1945   monoTriangulation2(upperVert, lowerVert, &vArray, 0, endIndex-beginIndex,
1946                      0, /*decreasing chain*/
1947                      pStream);
1948 }
1949
1950 /*return i, such that from begin to i-1 the chain is strictly u-monotone.
1951  */
1952 Int findIncreaseChainFromBegin(vertexArray* chain, Int begin ,Int end)
1953 {
1954   Int i=begin;
1955   Real prevU = chain->getVertex(i)[0];
1956   Real thisU;
1957   for(i=begin+1; i<=end; i++){
1958     thisU = chain->getVertex(i)[0];
1959
1960     if(prevU < thisU){
1961       prevU = thisU;
1962     }
1963     else
1964       break;
1965   }
1966   return i;
1967 }
1968
1969 /*check whether there is a vertex whose v value is strictly
1970  *inbetween vup vbelow
1971  *if no middle exists return -1, else return the idnex.
1972  */
1973 Int checkMiddle(vertexArray* chain, Int begin, Int end,
1974                 Real vup, Real vbelow)
1975 {
1976   Int i;
1977   for(i=begin; i<=end; i++)
1978     {
1979       if(chain->getVertex(i)[1] < vup && chain->getVertex(i)[1]>vbelow)
1980         return i;
1981     }
1982   return -1;
1983 }
1984
1985 /*the degenerat case of sampleLeftOneGridStep*/
1986 void sampleLeftOneGridStepNoMiddle(vertexArray* leftChain,
1987                                    Int beginLeftIndex,
1988                                    Int endLeftIndex,
1989                                    gridBoundaryChain* leftGridChain,
1990                                    Int leftGridChainStartIndex,
1991                                    primStream* pStream)
1992 {
1993   /*since there is no middle, there is at most one point which is on the
1994    *second grid line, there could be multiple points on the first (top)
1995    *grid line.
1996    */
1997
1998   leftGridChain->leftEndFan(leftGridChainStartIndex+1, pStream);
1999
2000   monoTriangulation2(leftGridChain->get_vertex(leftGridChainStartIndex),
2001                      leftGridChain->get_vertex(leftGridChainStartIndex+1),
2002                      leftChain,
2003                      beginLeftIndex,
2004                      endLeftIndex,
2005                      1, //is increase chain.
2006                      pStream);
2007 }
2008
2009
2010
2011 /*sampling the left area in between two grid lines.
2012  */
2013 void sampleLeftOneGridStep(vertexArray* leftChain,
2014                   Int beginLeftIndex,
2015                   Int endLeftIndex,
2016                   gridBoundaryChain* leftGridChain,
2017                   Int leftGridChainStartIndex,
2018                   primStream* pStream
2019                   )
2020 {
2021   if(checkMiddle(leftChain, beginLeftIndex, endLeftIndex,
2022                  leftGridChain->get_v_value(leftGridChainStartIndex),
2023                  leftGridChain->get_v_value(leftGridChainStartIndex+1))<0)
2024
2025     {
2026
2027       sampleLeftOneGridStepNoMiddle(leftChain, beginLeftIndex, endLeftIndex, leftGridChain, leftGridChainStartIndex, pStream);
2028       return;
2029     }
2030
2031   //copy into a polygon
2032   {
2033     directedLine* poly = NULL;
2034     sampledLine* sline;
2035     directedLine* dline;
2036     gridWrap* grid = leftGridChain->getGrid();
2037     Real vert1[2];
2038     Real vert2[2];
2039     Int i;
2040
2041     Int innerInd = leftGridChain->getInnerIndex(leftGridChainStartIndex+1);
2042     Int upperInd = leftGridChain->getUlineIndex(leftGridChainStartIndex);
2043     Int lowerInd = leftGridChain->getUlineIndex(leftGridChainStartIndex+1);
2044     Real upperV = leftGridChain->get_v_value(leftGridChainStartIndex);
2045     Real lowerV = leftGridChain->get_v_value(leftGridChainStartIndex+1);
2046
2047     //the upper gridline
2048     vert1[1] = vert2[1] = upperV;
2049     for(i=innerInd;     i>upperInd;     i--)
2050       {
2051         vert1[0]=grid->get_u_value(i);
2052         vert2[0]=grid->get_u_value(i-1);
2053         sline = new sampledLine(vert1, vert2);
2054         dline = new directedLine(INCREASING, sline);
2055         if(poly == NULL)
2056           poly = dline;
2057         else
2058           poly->insert(dline);
2059       }
2060
2061     //the edge connecting upper grid with left chain
2062     vert1[0] = grid->get_u_value(upperInd);
2063     vert1[1] = upperV;
2064     sline = new sampledLine(vert1, leftChain->getVertex(beginLeftIndex));
2065     dline = new directedLine(INCREASING, sline);
2066     if(poly == NULL)
2067       poly = dline;
2068     else
2069       poly->insert(dline);
2070
2071     //the left chain
2072     for(i=beginLeftIndex; i<endLeftIndex; i++)
2073       {
2074         sline = new sampledLine(leftChain->getVertex(i), leftChain->getVertex(i+1));
2075         dline = new directedLine(INCREASING, sline);
2076         poly->insert(dline);
2077       }
2078
2079     //the edge connecting left chain with lower gridline
2080     vert2[0] = grid->get_u_value(lowerInd);
2081     vert2[1] = lowerV;
2082     sline = new sampledLine(leftChain->getVertex(endLeftIndex), vert2);
2083     dline = new directedLine(INCREASING, sline);
2084     poly->insert(dline);
2085
2086     //the lower grid line
2087     vert1[1] = vert2[1] = lowerV;
2088     for(i=lowerInd; i<innerInd; i++)
2089       {
2090         vert1[0] = grid->get_u_value(i);
2091         vert2[0] = grid->get_u_value(i+1);
2092         sline = new sampledLine(vert1, vert2);
2093         dline = new directedLine(INCREASING, sline);
2094         poly->insert(dline);
2095       }
2096
2097     //the vertical grid line segement
2098     vert1[0]=vert2[0] = grid->get_u_value(innerInd);
2099     vert2[1]=upperV;
2100     vert1[1]=lowerV;
2101     sline=new sampledLine(vert1, vert2);
2102     dline=new directedLine(INCREASING, sline);
2103     poly->insert(dline);
2104     monoTriangulationOpt(poly, pStream);
2105     //cleanup
2106     poly->deleteSinglePolygonWithSline();
2107     return;
2108   }
2109
2110
2111
2112
2113
2114   Int i;
2115   if(1/*leftGridChain->getUlineIndex(leftGridChainStartIndex) >=
2116      leftGridChain->getUlineIndex(leftGridChainStartIndex+1)*/
2117      ) /*the second grid line is beyond the first one to the left*/
2118     {
2119       /*find the maximal U-monotone chain
2120        * of endLeftIndex, endLeftIndex-1, ...,
2121        */
2122       i=endLeftIndex;
2123       Real prevU = leftChain->getVertex(i)[0];
2124       for(i=endLeftIndex-1; i>=beginLeftIndex; i--){
2125         Real thisU = leftChain->getVertex(i)[0];
2126         if( prevU < thisU){
2127           prevU = thisU;
2128         }
2129         else
2130           break;
2131       }
2132       /*from endLeftIndex to i+1 is strictly U- monotone */
2133       /*if i+1==endLeftIndex and the vertex and leftchain is on the second gridline, then
2134        *we should use 2 vertices on the leftchain. If we only use one (endLeftIndex), then we
2135        *we would output degenerate triangles
2136        */
2137       if(i+1 == endLeftIndex && leftChain->getVertex(endLeftIndex)[1] == leftGridChain->get_v_value(1+leftGridChainStartIndex))
2138         i--;
2139
2140       Int j = beginLeftIndex/*endLeftIndex*/+1;
2141
2142
2143       if(leftGridChain->getInnerIndex(leftGridChainStartIndex+1) > leftGridChain->getUlineIndex(leftGridChainStartIndex))
2144         {
2145           j = findIncreaseChainFromBegin(leftChain, beginLeftIndex, i+1/*endLeftIndex*/);
2146
2147           Int temp = beginLeftIndex;
2148           /*now from begin to j-1 is strictly u-monotone*/
2149           /*if j-1 is on the first grid line, then we want to skip to the vertex which is strictly
2150            *below the grid line. This vertexmust exist since there is a 'corner turn' inbetween the two grid lines
2151            */
2152           if(j-1 == beginLeftIndex)
2153             {
2154               while(leftChain->getVertex(j-1)[1] == leftGridChain->get_v_value(leftGridChainStartIndex))
2155                 j++;
2156
2157               Real vert[2];
2158               vert[0] = leftGridChain->get_u_value(leftGridChainStartIndex);
2159               vert[1] = leftGridChain->get_v_value(leftGridChainStartIndex);
2160
2161               monoTriangulation2(
2162                                  vert/*leftChain->getVertex(beginLeftIndex)*/,
2163                                  leftChain->getVertex(j-1),
2164                                  leftChain,
2165                                  beginLeftIndex,
2166                                  j-2,
2167                                  1,
2168                                  pStream //increase chain
2169                                  );
2170
2171               temp = j-1;
2172             }
2173
2174           stripOfFanLeft(leftChain, j-1, temp/*beginLeftIndex*/, leftGridChain->getGrid(),
2175                          leftGridChain->getVlineIndex(leftGridChainStartIndex),
2176                          leftGridChain->getUlineIndex(leftGridChainStartIndex),
2177                          leftGridChain->getInnerIndex(leftGridChainStartIndex+1),
2178                          pStream,
2179                          1 /*the grid line is above the trim line*/
2180                          );
2181         }
2182
2183       stripOfFanLeft(leftChain, endLeftIndex, i+1, leftGridChain->getGrid(),
2184                      leftGridChain->getVlineIndex(leftGridChainStartIndex+1),
2185                      leftGridChain->getUlineIndex(leftGridChainStartIndex+1),
2186                      leftGridChain->getInnerIndex(leftGridChainStartIndex+1),
2187                      pStream,
2188                      0 /*the grid line is below the trim lines*/
2189                      );
2190
2191       /*monotone triangulate the remaining left chain togther with the
2192        *two vertices on the two grid v-lines.
2193        */
2194       Real vert[2][2];
2195       vert[0][0]=vert[1][0] = leftGridChain->getInner_u_value(leftGridChainStartIndex+1);
2196       vert[0][1] = leftGridChain->get_v_value(leftGridChainStartIndex);
2197       vert[1][1] = leftGridChain->get_v_value(leftGridChainStartIndex+1);
2198
2199 //      vertexArray right(vert, 2);
2200
2201       monoTriangulation2(
2202                          &vert[0][0], /*top vertex */
2203                          &vert[1][0], /*bottom vertex*/
2204                          leftChain,
2205                          /*beginLeftIndex*/j-1,
2206                          i+1,
2207                          1, /*an increasing chain*/
2208                          pStream);
2209     }
2210   else /*the second one is shorter than the first one to the left*/
2211     {
2212       /*find the maximal U-monotone chain of beginLeftIndex, beginLeftIndex+1,...,
2213        */
2214       i=beginLeftIndex;
2215       Real prevU = leftChain->getVertex(i)[0];
2216       for(i=beginLeftIndex+1; i<=endLeftIndex; i++){
2217         Real thisU = leftChain->getVertex(i)[0];
2218
2219         if(prevU < thisU){
2220           prevU = thisU;
2221         }
2222         else
2223           break;
2224       }
2225       /*from beginLeftIndex to i-1 is strictly U-monotone*/
2226
2227
2228       stripOfFanLeft(leftChain, i-1, beginLeftIndex, leftGridChain->getGrid(),
2229                          leftGridChain->getVlineIndex(leftGridChainStartIndex),
2230                          leftGridChain->getUlineIndex(leftGridChainStartIndex),
2231                          leftGridChain->getUlineIndex(leftGridChainStartIndex+1),
2232                          pStream,
2233                          1 /*the grid line is above the trim lines*/
2234                          );
2235       /*monotone triangulate the remaining left chain together with the
2236        *two vertices on the two grid v-lines.
2237        */
2238       Real vert[2][2];
2239       vert[0][0]=vert[1][0] = leftGridChain->get_u_value(leftGridChainStartIndex+1);
2240       vert[0][1] = leftGridChain->get_v_value(leftGridChainStartIndex);
2241       vert[1][1] = leftGridChain->get_v_value(leftGridChainStartIndex+1);
2242
2243       vertexArray right(vert, 2);
2244
2245       monoTriangulation2(
2246                          &vert[0][0], //top vertex
2247                          &vert[1][0], //bottom vertex
2248                          leftChain,
2249                          i-1,
2250                          endLeftIndex,
2251                          1, /*an increase chain*/
2252                          pStream);
2253
2254     }
2255 }
2256
2257 /*n_upper>=1
2258  *n_lower>=1
2259  */
2260 void triangulateXYMono(Int n_upper, Real upperVerts[][2],
2261                        Int n_lower, Real lowerVerts[][2],
2262                        primStream* pStream)
2263 {
2264   Int i,j,k,l;
2265   Real* leftMostV;
2266
2267   assert(n_upper>=1 && n_lower>=1);
2268   if(upperVerts[0][0] <= lowerVerts[0][0])
2269     {
2270       i=1;
2271       j=0;
2272       leftMostV = upperVerts[0];
2273     }
2274   else
2275     {
2276       i=0;
2277       j=1;
2278       leftMostV = lowerVerts[0];
2279     }
2280
2281   while(1)
2282     {
2283       if(i >= n_upper) /*case1: no more in upper*/
2284         {
2285
2286           if(j<n_lower-1) /*at least two vertices in lower*/
2287             {
2288               pStream->begin();
2289               pStream->insert(leftMostV);
2290               while(j<n_lower){
2291                 pStream->insert(lowerVerts[j]);
2292                 j++;
2293               }
2294               pStream->end(PRIMITIVE_STREAM_FAN);
2295             }
2296
2297           break;
2298         }
2299       else if(j>= n_lower) /*case2: no more in lower*/
2300         {
2301
2302           if(i<n_upper-1) /*at least two vertices in upper*/
2303             {
2304               pStream->begin();
2305               pStream->insert(leftMostV);
2306
2307               for(k=n_upper-1; k>=i; k--)
2308                 pStream->insert(upperVerts[k]);
2309
2310               pStream->end(PRIMITIVE_STREAM_FAN);
2311             }
2312
2313           break;
2314         }
2315       else /* case3: neither is empty, plus the leftMostV, there is at least one triangle to output*/
2316         {
2317
2318           if(upperVerts[i][0] <= lowerVerts[j][0])
2319             {
2320               pStream->begin();
2321               pStream->insert(lowerVerts[j]); /*the origin of this fan*/
2322
2323               /*find the last k>=i such that
2324                *upperverts[k][0] <= lowerverts[j][0]
2325                */
2326               k=i;
2327               while(k<n_upper)
2328                 {
2329                   if(upperVerts[k][0] > lowerVerts[j][0])
2330                     break;
2331                   k++;
2332                 }
2333               k--;
2334               for(l=k; l>=i; l--)/*the reverse is for two-face lighting*/
2335                 {
2336                   pStream->insert(upperVerts[l]);
2337                 }
2338               pStream->insert(leftMostV);
2339
2340               pStream->end(PRIMITIVE_STREAM_FAN);
2341               //update i for next loop
2342               i = k+1;
2343               leftMostV = upperVerts[k];
2344
2345             }
2346           else /*upperVerts[i][0] > lowerVerts[j][0]*/
2347             {
2348               pStream->begin();
2349               pStream->insert(upperVerts[i]);/*the origion of this fan*/
2350               pStream->insert(leftMostV);
2351               /*find the last k>=j such that
2352                *lowerverts[k][0] < upperverts[i][0]*/
2353               k=j;
2354               while(k< n_lower)
2355                 {
2356                   if(lowerVerts[k][0] >= upperVerts[i][0])
2357                     break;
2358                   pStream->insert(lowerVerts[k]);
2359                   k++;
2360                 }
2361               pStream->end(PRIMITIVE_STREAM_FAN);
2362               j=k;
2363               leftMostV = lowerVerts[j-1];
2364             }
2365         }
2366     }
2367 }
2368
2369
2370 void stripOfFanLeft(vertexArray* leftChain,
2371                     Int largeIndex,
2372                     Int smallIndex,
2373                     gridWrap* grid,
2374                     Int vlineIndex,
2375                     Int ulineSmallIndex,
2376                     Int ulineLargeIndex,
2377                     primStream* pStream,
2378                     Int gridLineUp /*1 if the grid line is above the trim lines*/
2379                      )
2380 {
2381   assert(largeIndex >= smallIndex);
2382
2383   Real grid_v_value;
2384   grid_v_value = grid->get_v_value(vlineIndex);
2385
2386   Real2* trimVerts=(Real2*) malloc(sizeof(Real2)* (largeIndex-smallIndex+1));
2387   assert(trimVerts);
2388
2389
2390   Real2* gridVerts=(Real2*) malloc(sizeof(Real2)* (ulineLargeIndex-ulineSmallIndex+1));
2391   assert(gridVerts);
2392
2393   Int k,i;
2394   if(gridLineUp) /*trim line is below grid line, so trim vertices are going right when index increases*/
2395     for(k=0, i=smallIndex; i<=largeIndex; i++, k++)
2396       {
2397       trimVerts[k][0] = leftChain->getVertex(i)[0];
2398       trimVerts[k][1] = leftChain->getVertex(i)[1];
2399     }
2400   else
2401     for(k=0, i=largeIndex; i>=smallIndex; i--, k++)
2402       {
2403         trimVerts[k][0] = leftChain->getVertex(i)[0];
2404         trimVerts[k][1] = leftChain->getVertex(i)[1];
2405       }
2406
2407   for(k=0, i=ulineSmallIndex; i<= ulineLargeIndex; i++, k++)
2408     {
2409       gridVerts[k][0] = grid->get_u_value(i);
2410       gridVerts[k][1] = grid_v_value;
2411     }
2412
2413   if(gridLineUp)
2414     triangulateXYMono(
2415                       ulineLargeIndex-ulineSmallIndex+1, gridVerts,
2416                       largeIndex-smallIndex+1, trimVerts,
2417                       pStream);
2418   else
2419     triangulateXYMono(largeIndex-smallIndex+1, trimVerts,
2420                       ulineLargeIndex-ulineSmallIndex+1, gridVerts,
2421                       pStream);
2422   free(trimVerts);
2423   free(gridVerts);
2424 }
2425
2426
2427
2428
2429