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