src/glu/sgi/libnurbs/internals/intersect.cc

   1 /*
   2 ** License Applicability. Except to the extent portions of this file are
   3 ** made subject to an alternative license as permitted in the SGI Free
   4 ** Software License B, Version 1.1 (the "License"), the contents of this
   5 ** file are subject only to the provisions of the License. You may not use
   6 ** this file except in compliance with the License. You may obtain a copy
   7 ** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
   8 ** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
   9 **
  10 ** http://oss.sgi.com/projects/FreeB
  11 **
  12 ** Note that, as provided in the License, the Software is distributed on an
  13 ** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
  14 ** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
  15 ** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
  16 ** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
  17 **
  18 ** Original Code. The Original Code is: OpenGL Sample Implementation,
  19 ** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
  20 ** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
  21 ** Copyright in any portions created by third parties is as indicated
  22 ** elsewhere herein. All Rights Reserved.
  23 **
  24 ** Additional Notice Provisions: The application programming interfaces
  25 ** established by SGI in conjunction with the Original Code are The
  26 ** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
  27 ** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
  28 ** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
  29 ** Window System(R) (Version 1.3), released October 19, 1998. This software
  30 ** was created using the OpenGL(R) version 1.2.1 Sample Implementation
  31 ** published by SGI, but has not been independently verified as being
  32 ** compliant with the OpenGL(R) version 1.2.1 Specification.
  33 */
  34
  35 /*
  36  * intersect.c++
  37  *
  38  * $Date: 2001/03/17 00:25:41 $ $Revision: 1.1 $
  39  * $Header: /home/krh/git/sync/mesa-cvs-repo/Mesa/src/glu/sgi/libnurbs/internals/intersect.cc,v 1.1 2001/03/17 00:25:41 brianp Exp $
  40  */
  41
  42 #include "glimports.h"
  43 #include "myassert.h"
  44 #include "mystdio.h"
  45 #include "subdivider.h"
  46 #include "arc.h"
  47 #include "bin.h"
  48 #include "backend.h"
  49 #include "trimvertpool.h"
  50
  51 /*#define NOTDEF*/
  52
  53 enum i_result { INTERSECT_VERTEX, INTERSECT_EDGE };
  54
  55 /* local functions */
  56 static int              arc_classify( Arc_ptr, int, REAL );
  57 static enum i_result    pwlarc_intersect( PwlArc *, int, REAL, int, int[3] );
  58
  59
  60 void
  61 Subdivider::partition( Bin & bin, Bin & left, Bin & intersections,
  62                 Bin & right, Bin & unknown, int param, REAL value )
  63 {
  64     Bin headonleft, headonright, tailonleft, tailonright;
  65
  66     for( Arc_ptr jarc = bin.removearc(); jarc; jarc = bin.removearc() ) {
  67
  68         REAL tdiff = jarc->tail()[param] - value;
  69         REAL hdiff = jarc->head()[param] - value;
  70
  71         if( tdiff > 0.0 ) {
  72             if( hdiff > 0.0 ) {
  73                 right.addarc( jarc  );
  74             } else if( hdiff == 0.0 ) {
  75                 tailonright.addarc( jarc  );
  76             } else {
  77                 Arc_ptr jtemp;
  78                 switch( arc_split(jarc, param, value, 0) ) {
  79                     case 2:
  80                         tailonright.addarc( jarc  );
  81                         headonleft.addarc( jarc->next  );
  82                         break;
  83                     case 31:
  84                         assert( jarc->head()[param] > value );
  85                         right.addarc( jarc  );
  86                         tailonright.addarc( jtemp = jarc->next  );
  87                         headonleft.addarc( jtemp->next  );
  88                         break;
  89                     case 32:
  90                         assert( jarc->head()[param] <= value );
  91                         tailonright .addarc( jarc  );
  92                         headonleft.addarc( jtemp = jarc->next  );
  93                         left.addarc( jtemp->next  );
  94                         break;
  95                     case 4:
  96                         right.addarc( jarc  );
  97                         tailonright.addarc( jtemp = jarc->next  );
  98                         headonleft.addarc( jtemp = jtemp->next  );
  99                         left.addarc( jtemp->next  );
 100                 }
 101             }
 102         } else if( tdiff == 0.0 ) {
 103             if( hdiff > 0.0 ) {
 104                 headonright.addarc( jarc  );
 105             } else if( hdiff == 0.0 ) {
 106                 unknown.addarc( jarc  );
 107             } else {
 108                 headonleft.addarc( jarc  );
 109             }
 110         } else {
 111             if( hdiff > 0.0 ) {
 112                 Arc_ptr jtemp;
 113                 switch( arc_split(jarc, param, value, 1) ) {
 114                     case 2:
 115                         tailonleft.addarc( jarc  );
 116                         headonright.addarc( jarc->next  );
 117                         break;
 118                     case 31:
 119                         assert( jarc->head()[param] < value );
 120                         left.addarc( jarc  );
 121                         tailonleft.addarc( jtemp = jarc->next  );
 122                         headonright.addarc( jtemp->next  );
 123                         break;
 124                     case 32:
 125                         assert( jarc->head()[param] >= value );
 126                         tailonleft.addarc( jarc  );
 127                         headonright.addarc( jtemp = jarc->next  );
 128                         right.addarc( jtemp->next  );
 129                         break;
 130                     case 4:
 131                         left.addarc( jarc  );
 132                         tailonleft.addarc( jtemp = jarc->next  );
 133                         headonright.addarc( jtemp = jtemp->next  );
 134                         right.addarc( jtemp->next  );
 135                 }
 136             } else if( hdiff == 0.0 ) {
 137                 tailonleft.addarc( jarc  );
 138             } else {
 139                 left.addarc( jarc  );
 140             }
 141         }
 142     }
 143     if( param == 0 ) {
 144         classify_headonleft_s( headonleft, intersections, left, value );
 145         classify_tailonleft_s( tailonleft, intersections, left, value );
 146         classify_headonright_s( headonright, intersections, right, value );
 147         classify_tailonright_s( tailonright, intersections, right, value );
 148     } else {
 149         classify_headonleft_t( headonleft, intersections, left, value );
 150         classify_tailonleft_t( tailonleft, intersections, left, value );
 151         classify_headonright_t( headonright, intersections, right, value );
 152         classify_tailonright_t( tailonright, intersections, right, value );
 153     }
 154 }
 155
 156 inline static void
 157 vert_interp( TrimVertex *n, TrimVertex *l, TrimVertex *r, int p, REAL val )
 158 {
 159     assert( val > l->param[p]);
 160     assert( val < r->param[p]);
 161
 162     n->nuid = l->nuid;
 163
 164     n->param[p] = val;
 165     if( l->param[1-p] != r->param[1-p]  ) {
 166         REAL ratio = (val - l->param[p]) / (r->param[p] - l->param[p]);
 167         n->param[1-p] = l->param[1-p] +
 168                         ratio * (r->param[1-p] - l->param[1-p]);
 169     } else {
 170         n->param[1-p] = l->param[1-p];
 171     }
 172 }
 173
 174 int
 175 Subdivider::arc_split( Arc_ptr jarc, int param, REAL value, int dir )
 176 {
 177     int         maxvertex = jarc->pwlArc->npts;
 178     Arc_ptr     jarc1, jarc2, jarc3;
 179     TrimVertex* v = jarc->pwlArc->pts;
 180
 181     int         loc[3];
 182     switch( pwlarc_intersect( jarc->pwlArc, param, value, dir, loc ) ) {
 183
 184                 // When the parameter value lands on a vertex, life is sweet
 185     case INTERSECT_VERTEX: {
 186             jarc1 = new(arcpool) Arc( jarc, new( pwlarcpool) PwlArc( maxvertex-loc[1], &v[loc[1]] ) );
 187             jarc->pwlArc->npts = loc[1] + 1;
 188             jarc1->next = jarc->next;
 189             jarc1->next->prev = jarc1;
 190             jarc->next = jarc1;
 191             jarc1->prev = jarc;
 192             assert(jarc->check() != 0);
 193             return 2;
 194         }
 195
 196                 // When the parameter value intersects an edge, we have to
 197                 // interpolate a new vertex.  There are special cases
 198                 // if the new vertex is adjacent to one or both of the
 199                 // endpoints of the arc.
 200     case INTERSECT_EDGE: {
 201             int i, j;
 202             if( dir == 0 ) {
 203                 i = loc[0];
 204                 j = loc[2];
 205             } else {
 206                 i = loc[2];
 207                 j = loc[0];
 208             }
 209
 210 #ifndef NOTDEF
 211             // The split is between vertices at index j and i, in that
 212             // order (j < i)
 213
 214             // JEB:  This code is my idea of how to do the split without
 215             // increasing the number of links.  I'm doing this so that
 216             // the is_rect routine can recognize rectangles created by
 217             // subdivision.  In exchange for simplifying the curve list,
 218             // however, it costs in allocated space and vertex copies.
 219
 220             TrimVertex *newjunk = trimvertexpool.get(maxvertex -i+1 /*-j*/);
 221             int k;
 222             for(k=0; k<maxvertex-i; k++)
 223               {
 224                 newjunk[k+1] = v[i+k];
 225                 newjunk[k+1].nuid = jarc->nuid;
 226               }
 227
 228             TrimVertex *vcopy = trimvertexpool.get(maxvertex);
 229             for(k=0; k<maxvertex; k++)
 230               {
 231                 vcopy[k].param[0] = v[k].param[0];
 232                 vcopy[k].param[1] = v[k].param[1];
 233               }
 234             jarc->pwlArc->pts=vcopy;
 235
 236             v[i].nuid = jarc->nuid;
 237             v[j].nuid = jarc->nuid;
 238             vert_interp( &newjunk[0], &v[loc[0]], &v[loc[2]], param, value );
 239
 240             if( showingDegenerate() )
 241                 backend.triangle( &v[i], &newjunk[0], &v[j] );
 242
 243             vcopy[j+1].param[0]=newjunk[0].param[0];
 244             vcopy[j+1].param[1]=newjunk[0].param[1];
 245
 246
 247             jarc1 = new(arcpool) Arc( jarc,
 248                         new(pwlarcpool) PwlArc(maxvertex-i+1 , newjunk ) );
 249
 250             jarc->pwlArc->npts = j+2;
 251             jarc1->next = jarc->next;
 252             jarc1->next->prev = jarc1;
 253             jarc->next = jarc1;
 254             jarc1->prev = jarc;
 255             assert(jarc->check() != 0);
 256
 257             return 2;
 258 #endif //not NOTDEF
 259                 // JEB: This is the original version:
 260 #ifdef NOTDEF
 261
 262             TrimVertex *newjunk = trimvertexpool.get(3);
 263             v[i].nuid = jarc->nuid;
 264             v[j].nuid = jarc->nuid;
 265             newjunk[0] = v[j];
 266             newjunk[2] = v[i];
 267             vert_interp( &newjunk[1], &v[loc[0]], &v[loc[2]], param, value );
 268
 269             if( showingDegenerate() )
 270                 backend.triangle( &newjunk[2], &newjunk[1], &newjunk[0] );
 271
 272                 // New vertex adjacent to both endpoints
 273             if (maxvertex == 2) {
 274                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 275                 jarc->pwlArc->npts = 2;
 276                 jarc->pwlArc->pts = newjunk;
 277                 jarc1->next = jarc->next;
 278                 jarc1->next->prev = jarc1;
 279                 jarc->next = jarc1;
 280                 jarc1->prev = jarc;
 281                 assert(jarc->check() != 0);
 282
 283                 return 2;
 284
 285                 // New vertex adjacent to ending point of arc
 286             } else if (maxvertex - j == 2) {
 287                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk ) );
 288                 jarc2 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 289                 jarc->pwlArc->npts = maxvertex-1;
 290                 jarc2->next = jarc->next;
 291                 jarc2->next->prev = jarc2;
 292                 jarc->next = jarc1;
 293                 jarc1->prev = jarc;
 294                 jarc1->next = jarc2;
 295                 jarc2->prev = jarc1;
 296                 assert(jarc->check() != 0);
 297                 return 31;
 298
 299                 // New vertex adjacent to starting point of arc
 300             } else if (i == 1) {
 301                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 302                 jarc2 = new(arcpool) Arc( jarc,
 303                         new(pwlarcpool) PwlArc( maxvertex-1, &jarc->pwlArc->pts[1] ) );
 304                 jarc->pwlArc->npts = 2;
 305                 jarc->pwlArc->pts = newjunk;
 306                 jarc2->next = jarc->next;
 307                 jarc2->next->prev = jarc2;
 308                 jarc->next = jarc1;
 309                 jarc1->prev = jarc;
 310                 jarc1->next = jarc2;
 311                 jarc2->prev = jarc1;
 312                 assert(jarc->check() != 0);
 313                 return 32;
 314
 315                 // It's somewhere in the middle
 316             } else {
 317                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk ) );
 318                 jarc2 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 319                 jarc3 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( maxvertex-i, v+i ) );
 320                 jarc->pwlArc->npts = j + 1;
 321                 jarc3->next = jarc->next;
 322                 jarc3->next->prev = jarc3;
 323                 jarc->next = jarc1;
 324                 jarc1->prev = jarc;
 325                 jarc1->next = jarc2;
 326                 jarc2->prev = jarc1;
 327                 jarc2->next = jarc3;
 328                 jarc3->prev = jarc2;
 329                 assert(jarc->check() != 0);
 330                 return 4;
 331             }
 332 #endif // NOTDEF
 333         }
 334         default:
 335         return -1; //picked -1 since it's not used
 336     }
 337 }
 338
 339 /*----------------------------------------------------------------------------
 340  * pwlarc_intersect -  find intersection of pwlArc and isoparametric line
 341  *----------------------------------------------------------------------------
 342  */
 343
 344 static enum i_result
 345 pwlarc_intersect(
 346     PwlArc *pwlArc,
 347     int param,
 348     REAL value,
 349     int dir,
 350     int loc[3] )
 351 {
 352     assert( pwlArc->npts > 0 );
 353
 354     if( dir ) {
 355         TrimVertex *v = pwlArc->pts;
 356         int imin = 0;
 357         int imax = pwlArc->npts - 1;
 358         assert( value > v[imin].param[param] );
 359         assert( value < v[imax].param[param] );
 360         while( (imax - imin) > 1 ) {
 361             int imid = (imax + imin)/2;
 362             if( v[imid].param[param] > value )
 363                 imax = imid;
 364             else if( v[imid].param[param] < value )
 365                 imin = imid;
 366             else {
 367                 loc[1] = imid;
 368                 return INTERSECT_VERTEX;
 369             }
 370         }
 371         loc[0] = imin;
 372         loc[2] = imax;
 373         return INTERSECT_EDGE;
 374     } else {
 375         TrimVertex *v = pwlArc->pts;
 376         int imax = 0;
 377         int imin = pwlArc->npts - 1;
 378         assert( value > v[imin].param[param] );
 379         assert( value < v[imax].param[param] );
 380         while( (imin - imax) > 1 ) {
 381             int imid = (imax + imin)/2;
 382             if( v[imid].param[param] > value )
 383                 imax = imid;
 384             else if( v[imid].param[param] < value )
 385                 imin = imid;
 386             else {
 387                 loc[1] = imid;
 388                 return INTERSECT_VERTEX;
 389             }
 390         }
 391         loc[0] = imin;
 392         loc[2] = imax;
 393         return INTERSECT_EDGE;
 394     }
 395 }
 396
 397 /*----------------------------------------------------------------------------
 398  * arc_classify - determine which side of a line a jarc lies
 399  *----------------------------------------------------------------------------
 400  */
 401
 402 static int
 403 arc_classify( Arc_ptr jarc, int param, REAL value )
 404 {
 405     REAL tdiff, hdiff;
 406     if( param == 0 ) {
 407         tdiff = jarc->tail()[0] - value;
 408         hdiff = jarc->head()[0] - value;
 409     } else {
 410         tdiff = jarc->tail()[1] - value;
 411         hdiff = jarc->head()[1] - value;
 412     }
 413
 414     if( tdiff > 0.0 ) {
 415         if( hdiff > 0.0 ) {
 416             return 0x11;
 417         } else if( hdiff == 0.0 ) {
 418             return 0x12;
 419         } else {
 420             return 0x10;
 421         }
 422     } else if( tdiff == 0.0 ) {
 423         if( hdiff > 0.0 ) {
 424             return 0x21;
 425         } else if( hdiff == 0.0 ) {
 426             return 0x22;
 427         } else {
 428             return 0x20;
 429         }
 430     } else {
 431         if( hdiff > 0.0 ) {
 432             return 0x01;
 433         } else if( hdiff == 0.0 ) {
 434             return 0x02;
 435         } else {
 436             return 0;
 437         }
 438     }
 439 }
 440
 441 void
 442 Subdivider::classify_tailonleft_s( Bin& bin, Bin& in, Bin& out, REAL val )
 443 {
 444     /* tail at left, head on line */
 445     Arc_ptr j;
 446
 447     while( (j = bin.removearc()) != NULL ) {
 448         assert( arc_classify( j, 0, val ) == 0x02 );
 449         j->clearitail();
 450
 451         REAL diff = j->next->head()[0] - val;
 452         if( diff > 0.0 ) {
 453             in.addarc( j );
 454         } else if( diff < 0.0 ) {
 455             if( ccwTurn_sl( j, j->next ) )
 456                 out.addarc( j );
 457             else
 458                 in.addarc( j );
 459         } else {
 460             if( j->next->tail()[1] > j->next->head()[1] )
 461                 in.addarc(j);
 462             else
 463                 out.addarc(j);
 464         }
 465     }
 466 }
 467
 468 void
 469 Subdivider::classify_tailonleft_t( Bin& bin, Bin& in, Bin& out, REAL val )
 470 {
 471     /* tail at left, head on line */
 472     Arc_ptr j;
 473
 474     while( (j = bin.removearc()) != NULL ) {
 475         assert( arc_classify( j, 1, val ) == 0x02 );
 476         j->clearitail();
 477
 478         REAL diff = j->next->head()[1] - val;
 479         if( diff > 0.0 ) {
 480             in.addarc( j );
 481         } else if( diff < 0.0 ) {
 482             if( ccwTurn_tl( j, j->next ) )
 483                 out.addarc( j );
 484             else
 485                 in.addarc( j );
 486         } else {
 487             if (j->next->tail()[0] > j->next->head()[0] )
 488                 out.addarc( j );
 489             else
 490                 in.addarc( j );
 491         }
 492     }
 493 }
 494
 495 void
 496 Subdivider::classify_headonleft_s( Bin& bin, Bin& in, Bin& out, REAL val )
 497 {
 498     /* tail on line, head at left */
 499     Arc_ptr j;
 500
 501     while( (j = bin.removearc()) != NULL ) {
 502         assert( arc_classify( j, 0, val ) == 0x20 );
 503
 504         j->setitail();
 505
 506         REAL diff = j->prev->tail()[0] - val;
 507         if( diff > 0.0 ) {
 508             out.addarc( j );
 509         } else if( diff < 0.0 ) {
 510             if( ccwTurn_sl( j->prev, j ) )
 511                 out.addarc( j );
 512             else
 513                 in.addarc( j );
 514         } else {
 515             if( j->prev->tail()[1] > j->prev->head()[1] )
 516                 in.addarc( j );
 517             else
 518                 out.addarc( j );
 519         }
 520     }
 521 }
 522
 523 void
 524 Subdivider::classify_headonleft_t( Bin& bin, Bin& in, Bin& out, REAL val )
 525 {
 526     /* tail on line, head at left */
 527     Arc_ptr j;
 528
 529     while( (j = bin.removearc()) != NULL ) {
 530         assert( arc_classify( j, 1, val ) == 0x20 );
 531         j->setitail();
 532
 533         REAL diff = j->prev->tail()[1] - val;
 534         if( diff > 0.0 ) {
 535             out.addarc( j );
 536         } else if( diff < 0.0 ) {
 537             if( ccwTurn_tl( j->prev, j ) )
 538                 out.addarc( j );
 539             else
 540                 in.addarc( j );
 541         } else {
 542             if( j->prev->tail()[0] > j->prev->head()[0] )
 543                 out.addarc( j );
 544             else
 545                 in.addarc( j );
 546         }
 547     }
 548 }
 549
 550
 551 void
 552 Subdivider::classify_tailonright_s( Bin& bin, Bin& in, Bin& out, REAL val )
 553 {
 554     /* tail at right, head on line */
 555     Arc_ptr j;
 556
 557     while( (j = bin.removearc()) != NULL ) {
 558         assert( arc_classify( j, 0, val ) == 0x12);
 559
 560         j->clearitail();
 561
 562         REAL diff = j->next->head()[0] - val;
 563         if( diff > 0.0 ) {
 564             if( ccwTurn_sr( j, j->next ) )
 565                 out.addarc( j );
 566             else
 567                 in.addarc( j );
 568         } else if( diff < 0.0 ) {
 569             in.addarc( j );
 570         } else {
 571             if( j->next->tail()[1] > j->next->head()[1] )
 572                 out.addarc( j );
 573             else
 574                 in.addarc( j );
 575         }
 576     }
 577 }
 578
 579 void
 580 Subdivider::classify_tailonright_t( Bin& bin, Bin& in, Bin& out, REAL val )
 581 {
 582     /* tail at right, head on line */
 583     Arc_ptr j;
 584
 585     while( (j = bin.removearc()) != NULL ) {
 586         assert( arc_classify( j, 1, val ) == 0x12);
 587
 588         j->clearitail();
 589
 590         REAL diff =  j->next->head()[1] - val;
 591         if( diff > 0.0 ) {
 592             if( ccwTurn_tr( j, j->next ) )
 593                 out.addarc( j );
 594             else
 595                 in.addarc( j );
 596         } else if( diff < 0.0 ) {
 597             in.addarc( j );
 598         } else {
 599             if( j->next->tail()[0] > j->next->head()[0] )
 600                 in.addarc( j );
 601             else
 602                 out.addarc( j );
 603         }
 604     }
 605 }
 606
 607 void
 608 Subdivider::classify_headonright_s( Bin& bin, Bin& in, Bin& out, REAL val )
 609 {
 610     /* tail on line, head at right */
 611     Arc_ptr j;
 612
 613     while( (j = bin.removearc()) != NULL ) {
 614         assert( arc_classify( j, 0, val ) == 0x21 );
 615
 616         j->setitail();
 617
 618         REAL diff = j->prev->tail()[0] - val;
 619         if( diff > 0.0 ) {
 620             if( ccwTurn_sr( j->prev, j ) )
 621                 out.addarc( j );
 622             else
 623                 in.addarc( j );
 624         } else if( diff < 0.0 ) {
 625             out.addarc( j );
 626         } else {
 627             if( j->prev->tail()[1] > j->prev->head()[1] )
 628                 out.addarc( j );
 629             else
 630                 in.addarc( j );
 631         }
 632     }
 633 }
 634
 635 void
 636 Subdivider::classify_headonright_t( Bin& bin, Bin& in, Bin& out, REAL val )
 637 {
 638     /* tail on line, head at right */
 639     Arc_ptr j;
 640
 641     while( (j = bin.removearc()) != NULL ) {
 642         assert( arc_classify( j, 1, val ) == 0x21 );
 643
 644         j->setitail();
 645
 646         REAL diff = j->prev->tail()[1] - val;
 647         if( diff > 0.0 ) {
 648             if( ccwTurn_tr( j->prev, j ) )
 649                 out.addarc( j );
 650             else
 651                 in.addarc( j );
 652         } else if( diff < 0.0 ) {
 653             out.addarc( j );
 654         } else {
 655             if( j->prev->tail()[0] > j->prev->head()[0] )
 656                 in.addarc( j );
 657             else
 658                 out.addarc( j );
 659         }
 660     }
 661 }
 662