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  */
  39
  40 #include "glimports.h"
  41 #include "myassert.h"
  42 #include "mystdio.h"
  43 #include "subdivider.h"
  44 #include "arc.h"
  45 #include "bin.h"
  46 #include "backend.h"
  47 #include "trimvertpool.h"
  48
  49 /*#define NOTDEF*/
  50
  51 enum i_result { INTERSECT_VERTEX, INTERSECT_EDGE };
  52
  53 /* local functions */
  54 #ifndef NDEBUG  // for asserts only
  55 static int              arc_classify( Arc_ptr, int, REAL );
  56 #endif
  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;
 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             Arc_ptr    jarc2, jarc3;
 262
 263             TrimVertex *newjunk = trimvertexpool.get(3);
 264             v[i].nuid = jarc->nuid;
 265             v[j].nuid = jarc->nuid;
 266             newjunk[0] = v[j];
 267             newjunk[2] = v[i];
 268             vert_interp( &newjunk[1], &v[loc[0]], &v[loc[2]], param, value );
 269
 270             if( showingDegenerate() )
 271                 backend.triangle( &newjunk[2], &newjunk[1], &newjunk[0] );
 272
 273                 // New vertex adjacent to both endpoints
 274             if (maxvertex == 2) {
 275                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 276                 jarc->pwlArc->npts = 2;
 277                 jarc->pwlArc->pts = newjunk;
 278                 jarc1->next = jarc->next;
 279                 jarc1->next->prev = jarc1;
 280                 jarc->next = jarc1;
 281                 jarc1->prev = jarc;
 282                 assert(jarc->check() != 0);
 283
 284                 return 2;
 285
 286                 // New vertex adjacent to ending point of arc
 287             } else if (maxvertex - j == 2) {
 288                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk ) );
 289                 jarc2 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 290                 jarc->pwlArc->npts = maxvertex-1;
 291                 jarc2->next = jarc->next;
 292                 jarc2->next->prev = jarc2;
 293                 jarc->next = jarc1;
 294                 jarc1->prev = jarc;
 295                 jarc1->next = jarc2;
 296                 jarc2->prev = jarc1;
 297                 assert(jarc->check() != 0);
 298                 return 31;
 299
 300                 // New vertex adjacent to starting point of arc
 301             } else if (i == 1) {
 302                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 303                 jarc2 = new(arcpool) Arc( jarc,
 304                         new(pwlarcpool) PwlArc( maxvertex-1, &jarc->pwlArc->pts[1] ) );
 305                 jarc->pwlArc->npts = 2;
 306                 jarc->pwlArc->pts = newjunk;
 307                 jarc2->next = jarc->next;
 308                 jarc2->next->prev = jarc2;
 309                 jarc->next = jarc1;
 310                 jarc1->prev = jarc;
 311                 jarc1->next = jarc2;
 312                 jarc2->prev = jarc1;
 313                 assert(jarc->check() != 0);
 314                 return 32;
 315
 316                 // It's somewhere in the middle
 317             } else {
 318                 jarc1 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk ) );
 319                 jarc2 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( 2, newjunk+1 ) );
 320                 jarc3 = new(arcpool) Arc( jarc, new(pwlarcpool) PwlArc( maxvertex-i, v+i ) );
 321                 jarc->pwlArc->npts = j + 1;
 322                 jarc3->next = jarc->next;
 323                 jarc3->next->prev = jarc3;
 324                 jarc->next = jarc1;
 325                 jarc1->prev = jarc;
 326                 jarc1->next = jarc2;
 327                 jarc2->prev = jarc1;
 328                 jarc2->next = jarc3;
 329                 jarc3->prev = jarc2;
 330                 assert(jarc->check() != 0);
 331                 return 4;
 332             }
 333 #endif // NOTDEF
 334         }
 335         default:
 336         return -1; //picked -1 since it's not used
 337     }
 338 }
 339
 340 /*----------------------------------------------------------------------------
 341  * pwlarc_intersect -  find intersection of pwlArc and isoparametric line
 342  *----------------------------------------------------------------------------
 343  */
 344
 345 static enum i_result
 346 pwlarc_intersect(
 347     PwlArc *pwlArc,
 348     int param,
 349     REAL value,
 350     int dir,
 351     int loc[3] )
 352 {
 353     assert( pwlArc->npts > 0 );
 354
 355     if( dir ) {
 356         TrimVertex *v = pwlArc->pts;
 357         int imin = 0;
 358         int imax = pwlArc->npts - 1;
 359         assert( value > v[imin].param[param] );
 360         assert( value < v[imax].param[param] );
 361         while( (imax - imin) > 1 ) {
 362             int imid = (imax + imin)/2;
 363             if( v[imid].param[param] > value )
 364                 imax = imid;
 365             else if( v[imid].param[param] < value )
 366                 imin = imid;
 367             else {
 368                 loc[1] = imid;
 369                 return INTERSECT_VERTEX;
 370             }
 371         }
 372         loc[0] = imin;
 373         loc[2] = imax;
 374         return INTERSECT_EDGE;
 375     } else {
 376         TrimVertex *v = pwlArc->pts;
 377         int imax = 0;
 378         int imin = pwlArc->npts - 1;
 379         assert( value > v[imin].param[param] );
 380         assert( value < v[imax].param[param] );
 381         while( (imin - imax) > 1 ) {
 382             int imid = (imax + imin)/2;
 383             if( v[imid].param[param] > value )
 384                 imax = imid;
 385             else if( v[imid].param[param] < value )
 386                 imin = imid;
 387             else {
 388                 loc[1] = imid;
 389                 return INTERSECT_VERTEX;
 390             }
 391         }
 392         loc[0] = imin;
 393         loc[2] = imax;
 394         return INTERSECT_EDGE;
 395     }
 396 }
 397
 398 /*----------------------------------------------------------------------------
 399  * arc_classify - determine which side of a line a jarc lies
 400  *----------------------------------------------------------------------------
 401  */
 402
 403 #ifndef NDEBUG  // for asserts only
 404 static int
 405 arc_classify( Arc_ptr jarc, int param, REAL value )
 406 {
 407     REAL tdiff, hdiff;
 408     if( param == 0 ) {
 409         tdiff = jarc->tail()[0] - value;
 410         hdiff = jarc->head()[0] - value;
 411     } else {
 412         tdiff = jarc->tail()[1] - value;
 413         hdiff = jarc->head()[1] - value;
 414     }
 415
 416     if( tdiff > 0.0 ) {
 417         if( hdiff > 0.0 ) {
 418             return 0x11;
 419         } else if( hdiff == 0.0 ) {
 420             return 0x12;
 421         } else {
 422             return 0x10;
 423         }
 424     } else if( tdiff == 0.0 ) {
 425         if( hdiff > 0.0 ) {
 426             return 0x21;
 427         } else if( hdiff == 0.0 ) {
 428             return 0x22;
 429         } else {
 430             return 0x20;
 431         }
 432     } else {
 433         if( hdiff > 0.0 ) {
 434             return 0x01;
 435         } else if( hdiff == 0.0 ) {
 436             return 0x02;
 437         } else {
 438             return 0;
 439         }
 440     }
 441 }
 442 #endif
 443
 444 void
 445 Subdivider::classify_tailonleft_s( Bin& bin, Bin& in, Bin& out, REAL val )
 446 {
 447     /* tail at left, head on line */
 448     Arc_ptr j;
 449
 450     while( (j = bin.removearc()) != NULL ) {
 451         assert( arc_classify( j, 0, val ) == 0x02 );
 452         j->clearitail();
 453
 454         REAL diff = j->next->head()[0] - val;
 455         if( diff > 0.0 ) {
 456             in.addarc( j );
 457         } else if( diff < 0.0 ) {
 458             if( ccwTurn_sl( j, j->next ) )
 459                 out.addarc( j );
 460             else
 461                 in.addarc( j );
 462         } else {
 463             if( j->next->tail()[1] > j->next->head()[1] )
 464                 in.addarc(j);
 465             else
 466                 out.addarc(j);
 467         }
 468     }
 469 }
 470
 471 void
 472 Subdivider::classify_tailonleft_t( Bin& bin, Bin& in, Bin& out, REAL val )
 473 {
 474     /* tail at left, head on line */
 475     Arc_ptr j;
 476
 477     while( (j = bin.removearc()) != NULL ) {
 478         assert( arc_classify( j, 1, val ) == 0x02 );
 479         j->clearitail();
 480
 481         REAL diff = j->next->head()[1] - val;
 482         if( diff > 0.0 ) {
 483             in.addarc( j );
 484         } else if( diff < 0.0 ) {
 485             if( ccwTurn_tl( j, j->next ) )
 486                 out.addarc( j );
 487             else
 488                 in.addarc( j );
 489         } else {
 490             if (j->next->tail()[0] > j->next->head()[0] )
 491                 out.addarc( j );
 492             else
 493                 in.addarc( j );
 494         }
 495     }
 496 }
 497
 498 void
 499 Subdivider::classify_headonleft_s( Bin& bin, Bin& in, Bin& out, REAL val )
 500 {
 501     /* tail on line, head at left */
 502     Arc_ptr j;
 503
 504     while( (j = bin.removearc()) != NULL ) {
 505         assert( arc_classify( j, 0, val ) == 0x20 );
 506
 507         j->setitail();
 508
 509         REAL diff = j->prev->tail()[0] - val;
 510         if( diff > 0.0 ) {
 511             out.addarc( j );
 512         } else if( diff < 0.0 ) {
 513             if( ccwTurn_sl( j->prev, j ) )
 514                 out.addarc( j );
 515             else
 516                 in.addarc( j );
 517         } else {
 518             if( j->prev->tail()[1] > j->prev->head()[1] )
 519                 in.addarc( j );
 520             else
 521                 out.addarc( j );
 522         }
 523     }
 524 }
 525
 526 void
 527 Subdivider::classify_headonleft_t( Bin& bin, Bin& in, Bin& out, REAL val )
 528 {
 529     /* tail on line, head at left */
 530     Arc_ptr j;
 531
 532     while( (j = bin.removearc()) != NULL ) {
 533         assert( arc_classify( j, 1, val ) == 0x20 );
 534         j->setitail();
 535
 536         REAL diff = j->prev->tail()[1] - val;
 537         if( diff > 0.0 ) {
 538             out.addarc( j );
 539         } else if( diff < 0.0 ) {
 540             if( ccwTurn_tl( j->prev, j ) )
 541                 out.addarc( j );
 542             else
 543                 in.addarc( j );
 544         } else {
 545             if( j->prev->tail()[0] > j->prev->head()[0] )
 546                 out.addarc( j );
 547             else
 548                 in.addarc( j );
 549         }
 550     }
 551 }
 552
 553
 554 void
 555 Subdivider::classify_tailonright_s( Bin& bin, Bin& in, Bin& out, REAL val )
 556 {
 557     /* tail at right, head on line */
 558     Arc_ptr j;
 559
 560     while( (j = bin.removearc()) != NULL ) {
 561         assert( arc_classify( j, 0, val ) == 0x12);
 562
 563         j->clearitail();
 564
 565         REAL diff = j->next->head()[0] - val;
 566         if( diff > 0.0 ) {
 567             if( ccwTurn_sr( j, j->next ) )
 568                 out.addarc( j );
 569             else
 570                 in.addarc( j );
 571         } else if( diff < 0.0 ) {
 572             in.addarc( j );
 573         } else {
 574             if( j->next->tail()[1] > j->next->head()[1] )
 575                 out.addarc( j );
 576             else
 577                 in.addarc( j );
 578         }
 579     }
 580 }
 581
 582 void
 583 Subdivider::classify_tailonright_t( Bin& bin, Bin& in, Bin& out, REAL val )
 584 {
 585     /* tail at right, head on line */
 586     Arc_ptr j;
 587
 588     while( (j = bin.removearc()) != NULL ) {
 589         assert( arc_classify( j, 1, val ) == 0x12);
 590
 591         j->clearitail();
 592
 593         REAL diff =  j->next->head()[1] - val;
 594         if( diff > 0.0 ) {
 595             if( ccwTurn_tr( j, j->next ) )
 596                 out.addarc( j );
 597             else
 598                 in.addarc( j );
 599         } else if( diff < 0.0 ) {
 600             in.addarc( j );
 601         } else {
 602             if( j->next->tail()[0] > j->next->head()[0] )
 603                 in.addarc( j );
 604             else
 605                 out.addarc( j );
 606         }
 607     }
 608 }
 609
 610 void
 611 Subdivider::classify_headonright_s( Bin& bin, Bin& in, Bin& out, REAL val )
 612 {
 613     /* tail on line, head at right */
 614     Arc_ptr j;
 615
 616     while( (j = bin.removearc()) != NULL ) {
 617         assert( arc_classify( j, 0, val ) == 0x21 );
 618
 619         j->setitail();
 620
 621         REAL diff = j->prev->tail()[0] - val;
 622         if( diff > 0.0 ) {
 623             if( ccwTurn_sr( j->prev, j ) )
 624                 out.addarc( j );
 625             else
 626                 in.addarc( j );
 627         } else if( diff < 0.0 ) {
 628             out.addarc( j );
 629         } else {
 630             if( j->prev->tail()[1] > j->prev->head()[1] )
 631                 out.addarc( j );
 632             else
 633                 in.addarc( j );
 634         }
 635     }
 636 }
 637
 638 void
 639 Subdivider::classify_headonright_t( Bin& bin, Bin& in, Bin& out, REAL val )
 640 {
 641     /* tail on line, head at right */
 642     Arc_ptr j;
 643
 644     while( (j = bin.removearc()) != NULL ) {
 645         assert( arc_classify( j, 1, val ) == 0x21 );
 646
 647         j->setitail();
 648
 649         REAL diff = j->prev->tail()[1] - val;
 650         if( diff > 0.0 ) {
 651             if( ccwTurn_tr( j->prev, j ) )
 652                 out.addarc( j );
 653             else
 654                 in.addarc( j );
 655         } else if( diff < 0.0 ) {
 656             out.addarc( j );
 657         } else {
 658             if( j->prev->tail()[0] > j->prev->head()[0] )
 659                 in.addarc( j );
 660             else
 661                 out.addarc( j );
 662         }
 663     }
 664 }
 665