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