lib/art/art_svp_wind.c

   1 /* Libart_LGPL - library of basic graphic primitives
   2  * Copyright (C) 1998-2000 Raph Levien
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Library General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Library General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Library General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19
  20 /* Primitive intersection and winding number operations on sorted
  21    vector paths.
  22
  23    These routines are internal to libart, used to construct operations
  24    like intersection, union, and difference. */
  25
  26 #include "config.h"
  27 #include "art_svp_wind.h"
  28
  29 #include <stdio.h> /* for printf of debugging info */
  30 #include <string.h> /* for memcpy */
  31 #include <math.h>
  32 #include "art_misc.h"
  33
  34 #include "art_rect.h"
  35 #include "art_svp.h"
  36
  37 #define noVERBOSE
  38
  39 #define PT_EQ(p1,p2) ((p1).x == (p2).x && (p1).y == (p2).y)
  40
  41 #define PT_CLOSE(p1,p2) (fabs ((p1).x - (p2).x) < 1e-6 && fabs ((p1).y - (p2).y) < 1e-6)
  42
  43 /* return nonzero and set *p to the intersection point if the lines
  44    z0-z1 and z2-z3 intersect each other. */
  45 static int
  46 intersect_lines (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3,
  47                  ArtPoint *p)
  48 {
  49   double a01, b01, c01;
  50   double a23, b23, c23;
  51   double d0, d1, d2, d3;
  52   double det;
  53
  54   /* if the vectors share an endpoint, they don't intersect */
  55   if (PT_EQ (z0, z2) || PT_EQ (z0, z3) || PT_EQ (z1, z2) || PT_EQ (z1, z3))
  56     return 0;
  57
  58 #if 0
  59   if (PT_CLOSE (z0, z2) || PT_CLOSE (z0, z3) || PT_CLOSE (z1, z2) || PT_CLOSE (z1, z3))
  60     return 0;
  61 #endif
  62
  63   /* find line equations ax + by + c = 0 */
  64   a01 = z0.y - z1.y;
  65   b01 = z1.x - z0.x;
  66   c01 = -(z0.x * a01 + z0.y * b01);
  67   /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y)
  68      = (z1.x * z0.y - z1.y * z0.x) */
  69
  70   d2 = a01 * z2.x + b01 * z2.y + c01;
  71   d3 = a01 * z3.x + b01 * z3.y + c01;
  72   if ((d2 > 0) == (d3 > 0))
  73     return 0;
  74
  75   a23 = z2.y - z3.y;
  76   b23 = z3.x - z2.x;
  77   c23 = -(z2.x * a23 + z2.y * b23);
  78
  79   d0 = a23 * z0.x + b23 * z0.y + c23;
  80   d1 = a23 * z1.x + b23 * z1.y + c23;
  81   if ((d0 > 0) == (d1 > 0))
  82     return 0;
  83
  84   /* now we definitely know that the lines intersect */
  85   /* solve the two linear equations ax + by + c = 0 */
  86   det = 1.0 / (a01 * b23 - a23 * b01);
  87   p->x = det * (c23 * b01 - c01 * b23);
  88   p->y = det * (c01 * a23 - c23 * a01);
  89
  90   return 1;
  91 }
  92
  93 #define EPSILON 1e-6
  94
  95 static double
  96 trap_epsilon (double v)
  97 {
  98   const double epsilon = EPSILON;
  99
 100   if (v < epsilon && v > -epsilon) return 0;
 101   else return v;
 102 }
 103
 104 /* Determine the order of line segments z0-z1 and z2-z3.
 105    Return +1 if z2-z3 lies entirely to the right of z0-z1,
 106    -1 if entirely to the left,
 107    or 0 if overlap.
 108
 109    The case analysis in this function is quite ugly. The fact that it's
 110    almost 200 lines long is ridiculous.
 111
 112    Ok, so here's the plan to cut it down:
 113
 114    First, do a bounding line comparison on the x coordinates. This is pretty
 115    much the common case, and should go quickly. It also takes care of the
 116    case where both lines are horizontal.
 117
 118    Then, do d0 and d1 computation, but only if a23 is nonzero.
 119
 120    Finally, do d2 and d3 computation, but only if a01 is nonzero.
 121
 122    Fall through to returning 0 (this will happen when both lines are
 123    horizontal and they overlap).
 124    */
 125 static int
 126 x_order (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3)
 127 {
 128   double a01, b01, c01;
 129   double a23, b23, c23;
 130   double d0, d1, d2, d3;
 131
 132   if (z0.y == z1.y)
 133     {
 134       if (z2.y == z3.y)
 135         {
 136           double x01min, x01max;
 137           double x23min, x23max;
 138
 139           if (z0.x > z1.x)
 140             {
 141               x01min = z1.x;
 142               x01max = z0.x;
 143             }
 144           else
 145             {
 146               x01min = z0.x;
 147               x01max = z1.x;
 148             }
 149
 150           if (z2.x > z3.x)
 151             {
 152               x23min = z3.x;
 153               x23max = z2.x;
 154             }
 155           else
 156             {
 157               x23min = z2.x;
 158               x23max = z3.x;
 159             }
 160
 161           if (x23min >= x01max) return 1;
 162           else if (x01min >= x23max) return -1;
 163           else return 0;
 164         }
 165       else
 166         {
 167           /* z0-z1 is horizontal, z2-z3 isn't */
 168           a23 = z2.y - z3.y;
 169           b23 = z3.x - z2.x;
 170           c23 = -(z2.x * a23 + z2.y * b23);
 171
 172           if (z3.y < z2.y)
 173             {
 174               a23 = -a23;
 175               b23 = -b23;
 176               c23 = -c23;
 177             }
 178
 179           d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23);
 180           d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23);
 181
 182           if (d0 > 0)
 183             {
 184               if (d1 >= 0) return 1;
 185               else return 0;
 186             }
 187           else if (d0 == 0)
 188             {
 189               if (d1 > 0) return 1;
 190               else if (d1 < 0) return -1;
 191               else /*printf ("case 1 degenerate\n")*/;
 192               return 0;
 193             }
 194           else /* d0 < 0 */
 195             {
 196               if (d1 <= 0) return -1;
 197               else return 0;
 198             }
 199         }
 200     }
 201   else if (z2.y == z3.y)
 202     {
 203       /* z2-z3 is horizontal, z0-z1 isn't */
 204       a01 = z0.y - z1.y;
 205       b01 = z1.x - z0.x;
 206       c01 = -(z0.x * a01 + z0.y * b01);
 207       /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y)
 208          = (z1.x * z0.y - z1.y * z0.x) */
 209
 210       if (z1.y < z0.y)
 211         {
 212           a01 = -a01;
 213           b01 = -b01;
 214           c01 = -c01;
 215         }
 216
 217       d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01);
 218       d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01);
 219
 220       if (d2 > 0)
 221         {
 222           if (d3 >= 0) return -1;
 223           else return 0;
 224         }
 225       else if (d2 == 0)
 226         {
 227           if (d3 > 0) return -1;
 228           else if (d3 < 0) return 1;
 229           else printf ("case 2 degenerate\n");
 230           return 0;
 231         }
 232       else /* d2 < 0 */
 233         {
 234           if (d3 <= 0) return 1;
 235           else return 0;
 236         }
 237     }
 238
 239   /* find line equations ax + by + c = 0 */
 240   a01 = z0.y - z1.y;
 241   b01 = z1.x - z0.x;
 242   c01 = -(z0.x * a01 + z0.y * b01);
 243   /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y)
 244      = -(z1.x * z0.y - z1.y * z0.x) */
 245
 246   if (a01 > 0)
 247     {
 248       a01 = -a01;
 249       b01 = -b01;
 250       c01 = -c01;
 251     }
 252   /* so now, (a01, b01) points to the left, thus a01 * x + b01 * y + c01
 253      is negative if the point lies to the right of the line */
 254
 255   d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01);
 256   d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01);
 257   if (d2 > 0)
 258     {
 259       if (d3 >= 0) return -1;
 260     }
 261   else if (d2 == 0)
 262     {
 263       if (d3 > 0) return -1;
 264       else if (d3 < 0) return 1;
 265       else
 266         /*fprintf (stderr, "colinear!\n")*/;
 267     }
 268   else /* d2 < 0 */
 269     {
 270       if (d3 <= 0) return 1;
 271     }
 272
 273   a23 = z2.y - z3.y;
 274   b23 = z3.x - z2.x;
 275   c23 = -(z2.x * a23 + z2.y * b23);
 276
 277   if (a23 > 0)
 278     {
 279       a23 = -a23;
 280       b23 = -b23;
 281       c23 = -c23;
 282     }
 283   d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23);
 284   d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23);
 285   if (d0 > 0)
 286     {
 287       if (d1 >= 0) return 1;
 288     }
 289   else if (d0 == 0)
 290     {
 291       if (d1 > 0) return 1;
 292       else if (d1 < 0) return -1;
 293       else
 294         /*fprintf (stderr, "colinear!\n")*/;
 295     }
 296   else /* d0 < 0 */
 297     {
 298       if (d1 <= 0) return -1;
 299     }
 300
 301   return 0;
 302 }
 303
 304 /* similar to x_order, but to determine whether point z0 + epsilon lies to
 305    the left of the line z2-z3 or to the right */
 306 static int
 307 x_order_2 (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3)
 308 {
 309   double a23, b23, c23;
 310   double d0, d1;
 311
 312   a23 = z2.y - z3.y;
 313   b23 = z3.x - z2.x;
 314   c23 = -(z2.x * a23 + z2.y * b23);
 315
 316   if (a23 > 0)
 317     {
 318       a23 = -a23;
 319       b23 = -b23;
 320       c23 = -c23;
 321     }
 322
 323   d0 = a23 * z0.x + b23 * z0.y + c23;
 324
 325   if (d0 > EPSILON)
 326     return -1;
 327   else if (d0 < -EPSILON)
 328     return 1;
 329
 330   d1 = a23 * z1.x + b23 * z1.y + c23;
 331   if (d1 > EPSILON)
 332     return -1;
 333   else if (d1 < -EPSILON)
 334     return 1;
 335
 336   if (z0.x == z1.x && z1.x == z2.x && z2.x == z3.x)
 337     {
 338       //art_dprint ("x_order_2: colinear and horizontally aligned!\n");
 339       return 0;
 340     }
 341
 342   if (z0.x <= z2.x && z1.x <= z2.x && z0.x <= z3.x && z1.x <= z3.x)
 343     return -1;
 344   if (z0.x >= z2.x && z1.x >= z2.x && z0.x >= z3.x && z1.x >= z3.x)
 345     return 1;
 346
 347   //fprintf (stderr, "x_order_2: colinear!\n");
 348   return 0;
 349 }
 350
 351 #ifdef DEAD_CODE
 352 /* Traverse the vector path, keeping it in x-sorted order.
 353
 354    This routine doesn't actually do anything - it's just here for
 355    explanatory purposes. */
 356 void
 357 traverse (ArtSVP *vp)
 358 {
 359   int *active_segs;
 360   int n_active_segs;
 361   int *cursor;
 362   int seg_idx;
 363   double y;
 364   int tmp1, tmp2;
 365   int asi;
 366   int i, j;
 367
 368   active_segs = art_new (int, vp->n_segs);
 369   cursor = art_new (int, vp->n_segs);
 370
 371   n_active_segs = 0;
 372   seg_idx = 0;
 373   y = vp->segs[0].points[0].y;
 374   while (seg_idx < vp->n_segs || n_active_segs > 0)
 375     {
 376       printf ("y = %g\n", y);
 377       /* delete segments ending at y from active list */
 378       for (i = 0; i < n_active_segs; i++)
 379         {
 380           asi = active_segs[i];
 381           if (vp->segs[asi].n_points - 1 == cursor[asi] &&
 382               vp->segs[asi].points[cursor[asi]].y == y)
 383             {
 384               printf ("deleting %d\n", asi);
 385               n_active_segs--;
 386               for (j = i; j < n_active_segs; j++)
 387                 active_segs[j] = active_segs[j + 1];
 388               i--;
 389             }
 390         }
 391
 392       /* insert new segments into the active list */
 393       while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
 394         {
 395           cursor[seg_idx] = 0;
 396           printf ("inserting %d\n", seg_idx);
 397           for (i = 0; i < n_active_segs; i++)
 398             {
 399               asi = active_segs[i];
 400               if (x_order (vp->segs[asi].points[cursor[asi]],
 401                            vp->segs[asi].points[cursor[asi] + 1],
 402                            vp->segs[seg_idx].points[0],
 403                            vp->segs[seg_idx].points[1]) == -1)
 404               break;
 405             }
 406           tmp1 = seg_idx;
 407           for (j = i; j < n_active_segs; j++)
 408             {
 409               tmp2 = active_segs[j];
 410               active_segs[j] = tmp1;
 411               tmp1 = tmp2;
 412             }
 413           active_segs[n_active_segs] = tmp1;
 414           n_active_segs++;
 415           seg_idx++;
 416         }
 417
 418       /* all active segs cross the y scanline (considering segs to be
 419        closed on top and open on bottom) */
 420       for (i = 0; i < n_active_segs; i++)
 421         {
 422           asi = active_segs[i];
 423           printf ("%d (%g, %g) - (%g, %g) %s\n", asi,
 424                   vp->segs[asi].points[cursor[asi]].x,
 425                   vp->segs[asi].points[cursor[asi]].y,
 426                   vp->segs[asi].points[cursor[asi] + 1].x,
 427                   vp->segs[asi].points[cursor[asi] + 1].y,
 428                   vp->segs[asi].dir ? "v" : "^");
 429         }
 430
 431       /* advance y to the next event */
 432       if (n_active_segs == 0)
 433         {
 434           if (seg_idx < vp->n_segs)
 435             y = vp->segs[seg_idx].points[0].y;
 436           /* else we're done */
 437         }
 438       else
 439         {
 440           asi = active_segs[0];
 441           y = vp->segs[asi].points[cursor[asi] + 1].y;
 442           for (i = 1; i < n_active_segs; i++)
 443             {
 444               asi = active_segs[i];
 445               if (y > vp->segs[asi].points[cursor[asi] + 1].y)
 446                 y = vp->segs[asi].points[cursor[asi] + 1].y;
 447             }
 448           if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
 449             y = vp->segs[seg_idx].points[0].y;
 450         }
 451
 452       /* advance cursors to reach new y */
 453       for (i = 0; i < n_active_segs; i++)
 454         {
 455           asi = active_segs[i];
 456           while (cursor[asi] < vp->segs[asi].n_points - 1 &&
 457                  y >= vp->segs[asi].points[cursor[asi] + 1].y)
 458             cursor[asi]++;
 459         }
 460       printf ("\n");
 461     }
 462   art_free (cursor);
 463   art_free (active_segs);
 464 }
 465 #endif
 466
 467 /* I believe that the loop will always break with i=1.
 468
 469    I think I'll want to change this from a simple sorted list to a
 470    modified stack. ips[*][0] will get its own data structure, and
 471    ips[*] will in general only be allocated if there is an intersection.
 472    Finally, the segment can be traced through the initial point
 473    (formerly ips[*][0]), backwards through the stack, and finally
 474    to cursor + 1.
 475
 476    This change should cut down on allocation bandwidth, and also
 477    eliminate the iteration through n_ipl below.
 478
 479 */
 480 static void
 481 insert_ip (int seg_i, int *n_ips, int *n_ips_max, ArtPoint **ips, ArtPoint ip)
 482 {
 483   int i;
 484   ArtPoint tmp1, tmp2;
 485   int n_ipl;
 486   ArtPoint *ipl;
 487
 488   n_ipl = n_ips[seg_i]++;
 489   if (n_ipl == n_ips_max[seg_i])
 490       art_expand (ips[seg_i], ArtPoint, n_ips_max[seg_i]);
 491   ipl = ips[seg_i];
 492   for (i = 1; i < n_ipl; i++)
 493     if (ipl[i].y > ip.y)
 494       break;
 495   tmp1 = ip;
 496   for (; i <= n_ipl; i++)
 497     {
 498       tmp2 = ipl[i];
 499       ipl[i] = tmp1;
 500       tmp1 = tmp2;
 501     }
 502 }
 503
 504 /* test active segment (i - 1) against i for intersection, if
 505    so, add intersection point to both ips lists. */
 506 static void
 507 intersect_neighbors (int i, int *active_segs,
 508                      int *n_ips, int *n_ips_max, ArtPoint **ips,
 509                      int *cursor, ArtSVP *vp)
 510 {
 511   ArtPoint z0, z1, z2, z3;
 512   int asi01, asi23;
 513   ArtPoint ip;
 514
 515   asi01 = active_segs[i - 1];
 516
 517   z0 = ips[asi01][0];
 518   if (n_ips[asi01] == 1)
 519     z1 = vp->segs[asi01].points[cursor[asi01] + 1];
 520   else
 521     z1 = ips[asi01][1];
 522
 523   asi23 = active_segs[i];
 524
 525   z2 = ips[asi23][0];
 526   if (n_ips[asi23] == 1)
 527     z3 = vp->segs[asi23].points[cursor[asi23] + 1];
 528   else
 529     z3 = ips[asi23][1];
 530
 531   if (intersect_lines (z0, z1, z2, z3, &ip))
 532     {
 533 #ifdef VERBOSE
 534       printf ("new intersection point: (%g, %g)\n", ip.x, ip.y);
 535 #endif
 536       insert_ip (asi01, n_ips, n_ips_max, ips, ip);
 537       insert_ip (asi23, n_ips, n_ips_max, ips, ip);
 538     }
 539 }
 540
 541 /* Add a new point to a segment in the svp.
 542
 543    Here, we also check to make sure that the segments satisfy nocross.
 544    However, this is only valuable for debugging, and could possibly be
 545    removed.
 546 */
 547 static void
 548 svp_add_point (ArtSVP *svp, int *n_points_max,
 549                ArtPoint p, int *seg_map, int *active_segs, int n_active_segs,
 550                int i)
 551 {
 552   int asi, asi_left, asi_right;
 553   int n_points, n_points_left, n_points_right;
 554   ArtSVPSeg *seg;
 555
 556   asi = seg_map[active_segs[i]];
 557   seg = &svp->segs[asi];
 558   n_points = seg->n_points;
 559   /* find out whether neighboring segments share a point */
 560   if (i > 0)
 561     {
 562       asi_left = seg_map[active_segs[i - 1]];
 563       n_points_left = svp->segs[asi_left].n_points;
 564       if (n_points_left > 1 &&
 565           PT_EQ (svp->segs[asi_left].points[n_points_left - 2],
 566                  svp->segs[asi].points[n_points - 1]))
 567         {
 568           /* ok, new vector shares a top point with segment to the left -
 569              now, check that it satisfies ordering invariant */
 570           if (x_order (svp->segs[asi_left].points[n_points_left - 2],
 571                        svp->segs[asi_left].points[n_points_left - 1],
 572                        svp->segs[asi].points[n_points - 1],
 573                        p) < 1)
 574
 575             {
 576 #ifdef VERBOSE
 577               printf ("svp_add_point: cross on left!\n");
 578 #endif
 579             }
 580         }
 581     }
 582
 583   if (i + 1 < n_active_segs)
 584     {
 585       asi_right = seg_map[active_segs[i + 1]];
 586       n_points_right = svp->segs[asi_right].n_points;
 587       if (n_points_right > 1 &&
 588           PT_EQ (svp->segs[asi_right].points[n_points_right - 2],
 589                  svp->segs[asi].points[n_points - 1]))
 590         {
 591           /* ok, new vector shares a top point with segment to the right -
 592              now, check that it satisfies ordering invariant */
 593           if (x_order (svp->segs[asi_right].points[n_points_right - 2],
 594                        svp->segs[asi_right].points[n_points_right - 1],
 595                        svp->segs[asi].points[n_points - 1],
 596                        p) > -1)
 597             {
 598 #ifdef VERBOSE
 599               printf ("svp_add_point: cross on right!\n");
 600 #endif
 601             }
 602         }
 603     }
 604   if (n_points_max[asi] == n_points)
 605     art_expand (seg->points, ArtPoint, n_points_max[asi]);
 606   seg->points[n_points] = p;
 607   if (p.x < seg->bbox.x0)
 608     seg->bbox.x0 = p.x;
 609   else if (p.x > seg->bbox.x1)
 610     seg->bbox.x1 = p.x;
 611   seg->bbox.y1 = p.y;
 612   seg->n_points++;
 613 }
 614
 615 #if 0
 616 /* find where the segment (currently at i) is supposed to go, and return
 617    the target index - if equal to i, then there is no crossing problem.
 618
 619    "Where it is supposed to go" is defined as following:
 620
 621    Delete element i, re-insert at position target (bumping everything
 622    target and greater to the right).
 623    */
 624 static int
 625 find_crossing (int i, int *active_segs, int n_active_segs,
 626                int *cursor, ArtPoint **ips, int *n_ips, ArtSVP *vp)
 627 {
 628   int asi, asi_left, asi_right;
 629   ArtPoint p0, p1;
 630   ArtPoint p0l, p1l;
 631   ArtPoint p0r, p1r;
 632   int target;
 633
 634   asi = active_segs[i];
 635   p0 = ips[asi][0];
 636   if (n_ips[asi] == 1)
 637     p1 = vp->segs[asi].points[cursor[asi] + 1];
 638   else
 639     p1 = ips[asi][1];
 640
 641   for (target = i; target > 0; target--)
 642     {
 643       asi_left = active_segs[target - 1];
 644       p0l = ips[asi_left][0];
 645       if (n_ips[asi_left] == 1)
 646         p1l = vp->segs[asi_left].points[cursor[asi_left] + 1];
 647       else
 648         p1l = ips[asi_left][1];
 649       if (!PT_EQ (p0, p0l))
 650         break;
 651
 652 #ifdef VERBOSE
 653       printf ("point matches on left (%g, %g) - (%g, %g) x (%g, %g) - (%g, %g)!\n",
 654               p0l.x, p0l.y, p1l.x, p1l.y, p0.x, p0.y, p1.x, p1.y);
 655 #endif
 656       if (x_order (p0l, p1l, p0, p1) == 1)
 657         break;
 658
 659 #ifdef VERBOSE
 660       printf ("scanning to the left (i=%d, target=%d)\n", i, target);
 661 #endif
 662     }
 663
 664   if (target < i) return target;
 665
 666   for (; target < n_active_segs - 1; target++)
 667     {
 668       asi_right = active_segs[target + 1];
 669       p0r = ips[asi_right][0];
 670       if (n_ips[asi_right] == 1)
 671         p1r = vp->segs[asi_right].points[cursor[asi_right] + 1];
 672       else
 673         p1r = ips[asi_right][1];
 674       if (!PT_EQ (p0, p0r))
 675         break;
 676
 677 #ifdef VERBOSE
 678       printf ("point matches on left (%g, %g) - (%g, %g) x (%g, %g) - (%g, %g)!\n",
 679               p0.x, p0.y, p1.x, p1.y, p0r.x, p0r.y, p1r.x, p1r.y);
 680 #endif
 681       if (x_order (p0r, p1r, p0, p1) == 1)
 682         break;
 683
 684 #ifdef VERBOSE
 685       printf ("scanning to the right (i=%d, target=%d)\n", i, target);
 686 #endif
 687     }
 688
 689   return target;
 690 }
 691 #endif
 692
 693 /* This routine handles the case where the segment changes its position
 694    in the active segment list. Generally, this will happen when the
 695    segment (defined by i and cursor) shares a top point with a neighbor,
 696    but breaks the ordering invariant.
 697
 698    Essentially, this routine sorts the lines [start..end), all of which
 699    share a top point. This is implemented as your basic insertion sort.
 700
 701    This routine takes care of intersecting the appropriate neighbors,
 702    as well.
 703
 704    A first argument of -1 immediately returns, which helps reduce special
 705    casing in the main unwind routine.
 706 */
 707 static void
 708 fix_crossing (int start, int end, int *active_segs, int n_active_segs,
 709               int *cursor, ArtPoint **ips, int *n_ips, int *n_ips_max,
 710               ArtSVP *vp, int *seg_map,
 711               ArtSVP **p_new_vp, int *pn_segs_max,
 712               int **pn_points_max)
 713 {
 714   int i, j;
 715   int target;
 716   int asi, asj;
 717   ArtPoint p0i, p1i;
 718   ArtPoint p0j, p1j;
 719   int swap = 0;
 720 #ifdef VERBOSE
 721   int k;
 722 #endif
 723   ArtPoint *pts;
 724
 725 #ifdef VERBOSE
 726   printf ("fix_crossing: [%d..%d)", start, end);
 727   for (k = 0; k < n_active_segs; k++)
 728     printf (" %d", active_segs[k]);
 729   printf ("\n");
 730 #endif
 731
 732   if (start == -1)
 733     return;
 734
 735   for (i = start + 1; i < end; i++)
 736     {
 737
 738       asi = active_segs[i];
 739       if (cursor[asi] < vp->segs[asi].n_points - 1) {
 740         p0i = ips[asi][0];
 741         if (n_ips[asi] == 1)
 742           p1i = vp->segs[asi].points[cursor[asi] + 1];
 743         else
 744           p1i = ips[asi][1];
 745
 746         for (j = i - 1; j >= start; j--)
 747           {
 748             asj = active_segs[j];
 749             if (cursor[asj] < vp->segs[asj].n_points - 1)
 750               {
 751                 p0j = ips[asj][0];
 752                 if (n_ips[asj] == 1)
 753                   p1j = vp->segs[asj].points[cursor[asj] + 1];
 754                 else
 755                   p1j = ips[asj][1];
 756
 757                 /* we _hope_ p0i = p0j */
 758                 if (x_order_2 (p0j, p1j, p0i, p1i) == -1)
 759                   break;
 760               }
 761           }
 762
 763         target = j + 1;
 764         /* target is where active_seg[i] _should_ be in active_segs */
 765
 766         if (target != i)
 767           {
 768             swap = 1;
 769
 770 #ifdef VERBOSE
 771             printf ("fix_crossing: at %i should be %i\n", i, target);
 772 #endif
 773
 774             /* let's close off all relevant segments */
 775             for (j = i; j >= target; j--)
 776               {
 777                 asi = active_segs[j];
 778                 /* First conjunct: this isn't the last point in the original
 779                    segment.
 780
 781                    Second conjunct: this isn't the first point in the new
 782                    segment (i.e. already broken).
 783                 */
 784                 if (cursor[asi] < vp->segs[asi].n_points - 1 &&
 785                     (*p_new_vp)->segs[seg_map[asi]].n_points != 1)
 786                   {
 787                     int seg_num;
 788                     /* so break here */
 789 #ifdef VERBOSE
 790                     printf ("closing off %d\n", j);
 791 #endif
 792
 793                     pts = art_new (ArtPoint, 16);
 794                     pts[0] = ips[asi][0];
 795                     seg_num = art_svp_add_segment (p_new_vp, pn_segs_max,
 796                                                    pn_points_max,
 797                                                    1, vp->segs[asi].dir,
 798                                                    pts,
 799                                                    NULL);
 800                     (*pn_points_max)[seg_num] = 16;
 801                     seg_map[asi] = seg_num;
 802                   }
 803               }
 804
 805             /* now fix the ordering in active_segs */
 806             asi = active_segs[i];
 807             for (j = i; j > target; j--)
 808               active_segs[j] = active_segs[j - 1];
 809             active_segs[j] = asi;
 810           }
 811       }
 812     }
 813   if (swap && start > 0)
 814     {
 815       int as_start;
 816
 817       as_start = active_segs[start];
 818       if (cursor[as_start] < vp->segs[as_start].n_points)
 819         {
 820 #ifdef VERBOSE
 821           printf ("checking intersection of %d, %d\n", start - 1, start);
 822 #endif
 823           intersect_neighbors (start, active_segs,
 824                                n_ips, n_ips_max, ips,
 825                                cursor, vp);
 826         }
 827     }
 828
 829   if (swap && end < n_active_segs)
 830     {
 831       int as_end;
 832
 833       as_end = active_segs[end - 1];
 834       if (cursor[as_end] < vp->segs[as_end].n_points)
 835         {
 836 #ifdef VERBOSE
 837           printf ("checking intersection of %d, %d\n", end - 1, end);
 838 #endif
 839           intersect_neighbors (end, active_segs,
 840                                n_ips, n_ips_max, ips,
 841                                cursor, vp);
 842         }
 843     }
 844   if (swap)
 845     {
 846 #ifdef VERBOSE
 847       printf ("fix_crossing return: [%d..%d)", start, end);
 848       for (k = 0; k < n_active_segs; k++)
 849         printf (" %d", active_segs[k]);
 850       printf ("\n");
 851 #endif
 852     }
 853 }
 854
 855 /* Return a new sorted vector that covers the same area as the
 856    argument, but which satisfies the nocross invariant.
 857
 858    Basically, this routine works by finding the intersection points,
 859    and cutting the segments at those points.
 860
 861    Status of this routine:
 862
 863    Basic correctness: Seems ok.
 864
 865    Numerical stability: known problems in the case of points falling
 866    on lines, and colinear lines. For actual use, randomly perturbing
 867    the vertices is currently recommended.
 868
 869    Speed: pretty good, although a more efficient priority queue, as
 870    well as bbox culling of potential intersections, are two
 871    optimizations that could help.
 872
 873    Precision: pretty good, although the numerical stability problems
 874    make this routine unsuitable for precise calculations of
 875    differences.
 876
 877 */
 878
 879 /* Here is a more detailed description of the algorithm. It follows
 880    roughly the structure of traverse (above), but is obviously quite
 881    a bit more complex.
 882
 883    Here are a few important data structures:
 884
 885    A new sorted vector path (new_svp).
 886
 887    For each (active) segment in the original, a list of intersection
 888    points.
 889
 890    Of course, the original being traversed.
 891
 892    The following invariants hold (in addition to the invariants
 893    of the traverse procedure).
 894
 895    The new sorted vector path lies entirely above the y scan line.
 896
 897    The new sorted vector path keeps the nocross invariant.
 898
 899    For each active segment, the y scan line crosses the line from the
 900    first to the second of the intersection points (where the second
 901    point is cursor + 1 if there is only one intersection point).
 902
 903    The list of intersection points + the (cursor + 1) point is kept
 904    in nondecreasing y order.
 905
 906    Of the active segments, none of the lines from first to second
 907    intersection point cross the 1st ip..2nd ip line of the left or
 908    right neighbor. (However, such a line may cross further
 909    intersection points of the neighbors, or segments past the
 910    immediate neighbors).
 911
 912    Of the active segments, all lines from 1st ip..2nd ip are in
 913    strictly increasing x_order (this is very similar to the invariant
 914    of the traverse procedure, but is explicitly stated here in terms
 915    of ips). (this basically says that nocross holds on the active
 916    segments)
 917
 918    The combination of the new sorted vector path, the path through all
 919    the intersection points to cursor + 1, and [cursor + 1, n_points)
 920    covers the same area as the argument.
 921
 922    Another important data structure is mapping from original segment
 923    number to new segment number.
 924
 925    The algorithm is perhaps best understood as advancing the cursors
 926    while maintaining these invariants. Here's roughly how it's done.
 927
 928    When deleting segments from the active list, those segments are added
 929    to the new sorted vector path. In addition, the neighbors may intersect
 930    each other, so they are intersection tested (see below).
 931
 932    When inserting new segments, they are intersection tested against
 933    their neighbors. The top point of the segment becomes the first
 934    intersection point.
 935
 936    Advancing the cursor is just a bit different from the traverse
 937    routine, as the cursor may advance through the intersection points
 938    as well. Only when there is a single intersection point in the list
 939    does the cursor advance in the original segment. In either case,
 940    the new vector is intersection tested against both neighbors. It
 941    also causes the vector over which the cursor is advancing to be
 942    added to the new svp.
 943
 944    Two steps need further clarification:
 945
 946    Intersection testing: the 1st ip..2nd ip lines of the neighbors
 947    are tested to see if they cross (using intersect_lines). If so,
 948    then the intersection point is added to the ip list of both
 949    segments, maintaining the invariant that the list of intersection
 950    points is nondecreasing in y).
 951
 952    Adding vector to new svp: if the new vector shares a top x
 953    coordinate with another vector, then it is checked to see whether
 954    it is in order. If not, then both segments are "broken," and then
 955    restarted. Note: in the case when both segments are in the same
 956    order, they may simply be swapped without breaking.
 957
 958    For the time being, I'm going to put some of these operations into
 959    subroutines. If it turns out to be a performance problem, I could
 960    try to reorganize the traverse procedure so that each is only
 961    called once, and inline them. But if it's not a performance
 962    problem, I'll just keep it this way, because it will probably help
 963    to make the code clearer, and I believe this code could use all the
 964    clarity it can get. */
 965 /**
 966  * art_svp_uncross: Resolve self-intersections of an svp.
 967  * @vp: The original svp.
 968  *
 969  * Finds all the intersections within @vp, and constructs a new svp
 970  * with new points added at these intersections.
 971  *
 972  * This routine needs to be redone from scratch with numerical robustness
 973  * in mind. I'm working on it.
 974  *
 975  * Return value: The new svp.
 976  **/
 977 ArtSVP *
 978 art_svp_uncross (ArtSVP *vp)
 979 {
 980   int *active_segs;
 981   int n_active_segs;
 982   int *cursor;
 983   int seg_idx;
 984   double y;
 985   int tmp1, tmp2;
 986   int asi;
 987   int i, j;
 988   /* new data structures */
 989   /* intersection points; invariant: *ips[i] is only allocated if
 990      i is active */
 991   int *n_ips, *n_ips_max;
 992   ArtPoint **ips;
 993   /* new sorted vector path */
 994   int n_segs_max, seg_num;
 995   ArtSVP *new_vp;
 996   int *n_points_max;
 997   /* mapping from argument to new segment numbers - again, only valid
 998    if active */
 999   int *seg_map;
1000   double y_curs;
1001   ArtPoint p_curs;
1002   int first_share;
1003   double share_x;
1004   ArtPoint *pts;
1005
1006   n_segs_max = 16;
1007   new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) +
1008                                 (n_segs_max - 1) * sizeof(ArtSVPSeg));
1009   new_vp->n_segs = 0;
1010
1011   if (vp->n_segs == 0)
1012     return new_vp;
1013
1014   active_segs = art_new (int, vp->n_segs);
1015   cursor = art_new (int, vp->n_segs);
1016
1017   seg_map = art_new (int, vp->n_segs);
1018   n_ips = art_new (int, vp->n_segs);
1019   n_ips_max = art_new (int, vp->n_segs);
1020   ips = art_new (ArtPoint *, vp->n_segs);
1021
1022   n_points_max = art_new (int, n_segs_max);
1023
1024   n_active_segs = 0;
1025   seg_idx = 0;
1026   y = vp->segs[0].points[0].y;
1027   while (seg_idx < vp->n_segs || n_active_segs > 0)
1028     {
1029 #ifdef VERBOSE
1030       printf ("y = %g\n", y);
1031 #endif
1032
1033       /* maybe move deletions to end of loop (to avoid so much special
1034          casing on the end of a segment)? */
1035
1036       /* delete segments ending at y from active list */
1037       for (i = 0; i < n_active_segs; i++)
1038         {
1039           asi = active_segs[i];
1040           if (vp->segs[asi].n_points - 1 == cursor[asi] &&
1041               vp->segs[asi].points[cursor[asi]].y == y)
1042             {
1043               do
1044                 {
1045 #ifdef VERBOSE
1046                   printf ("deleting %d\n", asi);
1047 #endif
1048                   art_free (ips[asi]);
1049                   n_active_segs--;
1050                   for (j = i; j < n_active_segs; j++)
1051                     active_segs[j] = active_segs[j + 1];
1052                   if (i < n_active_segs)
1053                     asi = active_segs[i];
1054                   else
1055                     break;
1056                 }
1057               while (vp->segs[asi].n_points - 1 == cursor[asi] &&
1058                      vp->segs[asi].points[cursor[asi]].y == y);
1059
1060               /* test intersection of neighbors */
1061               if (i > 0 && i < n_active_segs)
1062                 intersect_neighbors (i, active_segs,
1063                                      n_ips, n_ips_max, ips,
1064                                      cursor, vp);
1065
1066               i--;
1067             }
1068         }
1069
1070       /* insert new segments into the active list */
1071       while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
1072         {
1073 #ifdef VERBOSE
1074           printf ("inserting %d\n", seg_idx);
1075 #endif
1076           cursor[seg_idx] = 0;
1077           for (i = 0; i < n_active_segs; i++)
1078             {
1079               asi = active_segs[i];
1080               if (x_order_2 (vp->segs[seg_idx].points[0],
1081                              vp->segs[seg_idx].points[1],
1082                              vp->segs[asi].points[cursor[asi]],
1083                              vp->segs[asi].points[cursor[asi] + 1]) == -1)
1084                 break;
1085             }
1086
1087           /* Create and initialize the intersection points data structure */
1088           n_ips[seg_idx] = 1;
1089           n_ips_max[seg_idx] = 2;
1090           ips[seg_idx] = art_new (ArtPoint, n_ips_max[seg_idx]);
1091           ips[seg_idx][0] = vp->segs[seg_idx].points[0];
1092
1093           /* Start a new segment in the new vector path */
1094           pts = art_new (ArtPoint, 16);
1095           pts[0] = vp->segs[seg_idx].points[0];
1096           seg_num = art_svp_add_segment (&new_vp, &n_segs_max,
1097                                          &n_points_max,
1098                                          1, vp->segs[seg_idx].dir,
1099                                          pts,
1100                                          NULL);
1101           n_points_max[seg_num] = 16;
1102           seg_map[seg_idx] = seg_num;
1103
1104           tmp1 = seg_idx;
1105           for (j = i; j < n_active_segs; j++)
1106             {
1107               tmp2 = active_segs[j];
1108               active_segs[j] = tmp1;
1109               tmp1 = tmp2;
1110             }
1111           active_segs[n_active_segs] = tmp1;
1112           n_active_segs++;
1113
1114           if (i > 0)
1115             intersect_neighbors (i, active_segs,
1116                                  n_ips, n_ips_max, ips,
1117                                  cursor, vp);
1118
1119           if (i + 1 < n_active_segs)
1120             intersect_neighbors (i + 1, active_segs,
1121                                  n_ips, n_ips_max, ips,
1122                                  cursor, vp);
1123
1124           seg_idx++;
1125         }
1126
1127       /* all active segs cross the y scanline (considering segs to be
1128        closed on top and open on bottom) */
1129 #ifdef VERBOSE
1130       for (i = 0; i < n_active_segs; i++)
1131         {
1132           asi = active_segs[i];
1133           printf ("%d ", asi);
1134           for (j = 0; j < n_ips[asi]; j++)
1135             printf ("(%g, %g) - ",
1136                     ips[asi][j].x,
1137                     ips[asi][j].y);
1138           printf ("(%g, %g) %s\n",
1139                   vp->segs[asi].points[cursor[asi] + 1].x,
1140                   vp->segs[asi].points[cursor[asi] + 1].y,
1141                   vp->segs[asi].dir ? "v" : "^");
1142         }
1143 #endif
1144
1145       /* advance y to the next event
1146        Note: this is quadratic. We'd probably get decent constant
1147        factor speed improvement by caching the y_curs values. */
1148       if (n_active_segs == 0)
1149         {
1150           if (seg_idx < vp->n_segs)
1151             y = vp->segs[seg_idx].points[0].y;
1152           /* else we're done */
1153         }
1154       else
1155         {
1156           asi = active_segs[0];
1157           if (n_ips[asi] == 1)
1158             y = vp->segs[asi].points[cursor[asi] + 1].y;
1159           else
1160             y = ips[asi][1].y;
1161           for (i = 1; i < n_active_segs; i++)
1162             {
1163               asi = active_segs[i];
1164               if (n_ips[asi] == 1)
1165                 y_curs = vp->segs[asi].points[cursor[asi] + 1].y;
1166               else
1167                 y_curs = ips[asi][1].y;
1168               if (y > y_curs)
1169                 y = y_curs;
1170             }
1171           if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
1172             y = vp->segs[seg_idx].points[0].y;
1173         }
1174
1175       first_share = -1;
1176       share_x = 0; /* to avoid gcc warning, although share_x is never
1177                       used when first_share is -1 */
1178       /* advance cursors to reach new y */
1179       for (i = 0; i < n_active_segs; i++)
1180         {
1181           asi = active_segs[i];
1182           if (n_ips[asi] == 1)
1183             p_curs = vp->segs[asi].points[cursor[asi] + 1];
1184           else
1185             p_curs = ips[asi][1];
1186           if (p_curs.y == y)
1187             {
1188               svp_add_point (new_vp, n_points_max,
1189                              p_curs, seg_map, active_segs, n_active_segs, i);
1190
1191               n_ips[asi]--;
1192               for (j = 0; j < n_ips[asi]; j++)
1193                 ips[asi][j] = ips[asi][j + 1];
1194
1195               if (n_ips[asi] == 0)
1196                 {
1197                   ips[asi][0] = p_curs;
1198                   n_ips[asi] = 1;
1199                   cursor[asi]++;
1200                 }
1201
1202               if (first_share < 0 || p_curs.x != share_x)
1203                 {
1204                   /* this is where crossings are detected, and if
1205                      found, the active segments switched around. */
1206
1207                   fix_crossing (first_share, i,
1208                                 active_segs, n_active_segs,
1209                                 cursor, ips, n_ips, n_ips_max, vp, seg_map,
1210                                 &new_vp,
1211                                 &n_segs_max, &n_points_max);
1212
1213                   first_share = i;
1214                   share_x = p_curs.x;
1215                 }
1216
1217               if (cursor[asi] < vp->segs[asi].n_points - 1)
1218                 {
1219
1220                   if (i > 0)
1221                     intersect_neighbors (i, active_segs,
1222                                          n_ips, n_ips_max, ips,
1223                                          cursor, vp);
1224
1225                   if (i + 1 < n_active_segs)
1226                     intersect_neighbors (i + 1, active_segs,
1227                                          n_ips, n_ips_max, ips,
1228                                          cursor, vp);
1229                 }
1230             }
1231           else
1232             {
1233               /* not on a cursor point */
1234               fix_crossing (first_share, i,
1235                             active_segs, n_active_segs,
1236                             cursor, ips, n_ips, n_ips_max, vp, seg_map,
1237                             &new_vp,
1238                             &n_segs_max, &n_points_max);
1239               first_share = -1;
1240             }
1241         }
1242
1243       /* fix crossing on last shared group */
1244       fix_crossing (first_share, i,
1245                     active_segs, n_active_segs,
1246                     cursor, ips, n_ips, n_ips_max, vp, seg_map,
1247                     &new_vp,
1248                     &n_segs_max, &n_points_max);
1249
1250 #ifdef VERBOSE
1251       printf ("\n");
1252 #endif
1253     }
1254
1255   /* not necessary to sort, new segments only get added at y, which
1256      increases monotonically */
1257 #if 0
1258   qsort (&new_vp->segs, new_vp->n_segs, sizeof (svp_seg), svp_seg_compare);
1259   {
1260     int k;
1261     for (k = 0; k < new_vp->n_segs - 1; k++)
1262       {
1263         printf ("(%g, %g) - (%g, %g) %s (%g, %g) - (%g, %g)\n",
1264                 new_vp->segs[k].points[0].x,
1265                 new_vp->segs[k].points[0].y,
1266                 new_vp->segs[k].points[1].x,
1267                 new_vp->segs[k].points[1].y,
1268                 svp_seg_compare (&new_vp->segs[k], &new_vp->segs[k + 1]) > 1 ? ">": "<",
1269                 new_vp->segs[k + 1].points[0].x,
1270                 new_vp->segs[k + 1].points[0].y,
1271                 new_vp->segs[k + 1].points[1].x,
1272                 new_vp->segs[k + 1].points[1].y);
1273       }
1274   }
1275 #endif
1276
1277   art_free (n_points_max);
1278   art_free (seg_map);
1279   art_free (n_ips_max);
1280   art_free (n_ips);
1281   art_free (ips);
1282   art_free (cursor);
1283   art_free (active_segs);
1284
1285   return new_vp;
1286 }
1287
1288 #define noVERBOSE
1289
1290 /* Rewind a svp satisfying the nocross invariant.
1291
1292    The winding number of a segment is defined as the winding number of
1293    the points to the left while travelling in the direction of the
1294    segment. Therefore it preincrements and postdecrements as a scan
1295    line is traversed from left to right.
1296
1297    Status of this routine:
1298
1299    Basic correctness: Was ok in gfonted. However, this code does not
1300    yet compute bboxes for the resulting svp segs.
1301
1302    Numerical stability: known problems in the case of horizontal
1303    segments in polygons with any complexity. For actual use, randomly
1304    perturbing the vertices is recommended.
1305
1306    Speed: good.
1307
1308    Precision: good, except that no attempt is made to remove "hair".
1309    Doing random perturbation just makes matters worse.
1310
1311 */
1312 /**
1313  * art_svp_rewind_uncrossed: Rewind an svp satisfying the nocross invariant.
1314  * @vp: The original svp.
1315  * @rule: The winding rule.
1316  *
1317  * Creates a new svp with winding number of 0 or 1 everywhere. The @rule
1318  * argument specifies a rule for how winding numbers in the original
1319  * @vp map to the winding numbers in the result.
1320  *
1321  * With @rule == ART_WIND_RULE_NONZERO, the resulting svp has a
1322  * winding number of 1 where @vp has a nonzero winding number.
1323  *
1324  * With @rule == ART_WIND_RULE_INTERSECT, the resulting svp has a
1325  * winding number of 1 where @vp has a winding number greater than
1326  * 1. It is useful for computing intersections.
1327  *
1328  * With @rule == ART_WIND_RULE_ODDEVEN, the resulting svp has a
1329  * winding number of 1 where @vp has an odd winding number. It is
1330  * useful for implementing the even-odd winding rule of the
1331  * PostScript imaging model.
1332  *
1333  * With @rule == ART_WIND_RULE_POSITIVE, the resulting svp has a
1334  * winding number of 1 where @vp has a positive winding number. It is
1335  * useful for implementing asymmetric difference.
1336  *
1337  * This routine needs to be redone from scratch with numerical robustness
1338  * in mind. I'm working on it.
1339  *
1340  * Return value: The new svp.
1341  **/
1342 ArtSVP *
1343 art_svp_rewind_uncrossed (ArtSVP *vp, ArtWindRule rule)
1344 {
1345   int *active_segs;
1346   int n_active_segs;
1347   int *cursor;
1348   int seg_idx;
1349   double y;
1350   int tmp1, tmp2;
1351   int asi;
1352   int i, j;
1353
1354   ArtSVP *new_vp;
1355   int n_segs_max;
1356   int *winding;
1357   int left_wind;
1358   int wind;
1359   int keep, invert;
1360
1361 #ifdef VERBOSE
1362   print_svp (vp);
1363 #endif
1364   n_segs_max = 16;
1365   new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) +
1366                                 (n_segs_max - 1) * sizeof(ArtSVPSeg));
1367   new_vp->n_segs = 0;
1368
1369   if (vp->n_segs == 0)
1370     return new_vp;
1371
1372   winding = art_new (int, vp->n_segs);
1373
1374   active_segs = art_new (int, vp->n_segs);
1375   cursor = art_new (int, vp->n_segs);
1376
1377   n_active_segs = 0;
1378   seg_idx = 0;
1379   y = vp->segs[0].points[0].y;
1380   while (seg_idx < vp->n_segs || n_active_segs > 0)
1381     {
1382 #ifdef VERBOSE
1383       printf ("y = %g\n", y);
1384 #endif
1385       /* delete segments ending at y from active list */
1386       for (i = 0; i < n_active_segs; i++)
1387         {
1388           asi = active_segs[i];
1389           if (vp->segs[asi].n_points - 1 == cursor[asi] &&
1390               vp->segs[asi].points[cursor[asi]].y == y)
1391             {
1392 #ifdef VERBOSE
1393               printf ("deleting %d\n", asi);
1394 #endif
1395               n_active_segs--;
1396               for (j = i; j < n_active_segs; j++)
1397                 active_segs[j] = active_segs[j + 1];
1398               i--;
1399             }
1400         }
1401
1402       /* insert new segments into the active list */
1403       while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
1404         {
1405 #ifdef VERBOSE
1406           printf ("inserting %d\n", seg_idx);
1407 #endif
1408           cursor[seg_idx] = 0;
1409           for (i = 0; i < n_active_segs; i++)
1410             {
1411               asi = active_segs[i];
1412               if (x_order_2 (vp->segs[seg_idx].points[0],
1413                              vp->segs[seg_idx].points[1],
1414                              vp->segs[asi].points[cursor[asi]],
1415                              vp->segs[asi].points[cursor[asi] + 1]) == -1)
1416                 break;
1417             }
1418
1419           /* Determine winding number for this segment */
1420           if (i == 0)
1421             left_wind = 0;
1422           else if (vp->segs[active_segs[i - 1]].dir)
1423             left_wind = winding[active_segs[i - 1]];
1424           else
1425             left_wind = winding[active_segs[i - 1]] - 1;
1426
1427           if (vp->segs[seg_idx].dir)
1428             wind = left_wind + 1;
1429           else
1430             wind = left_wind;
1431
1432           winding[seg_idx] = wind;
1433
1434           switch (rule)
1435             {
1436             case ART_WIND_RULE_NONZERO:
1437               keep = (wind == 1 || wind == 0);
1438               invert = (wind == 0);
1439               break;
1440             case ART_WIND_RULE_INTERSECT:
1441               keep = (wind == 2);
1442               invert = 0;
1443               break;
1444             case ART_WIND_RULE_ODDEVEN:
1445               keep = 1;
1446               invert = !(wind & 1);
1447               break;
1448             case ART_WIND_RULE_POSITIVE:
1449               keep = (wind == 1);
1450               invert = 0;
1451               break;
1452             default:
1453               keep = 0;
1454               invert = 0;
1455               break;
1456             }
1457
1458           if (keep)
1459             {
1460               ArtPoint *points, *new_points;
1461               int n_points;
1462               int new_dir;
1463
1464 #ifdef VERBOSE
1465               printf ("keeping segment %d\n", seg_idx);
1466 #endif
1467               n_points = vp->segs[seg_idx].n_points;
1468               points = vp->segs[seg_idx].points;
1469               new_points = art_new (ArtPoint, n_points);
1470               memcpy (new_points, points, n_points * sizeof (ArtPoint));
1471               new_dir = vp->segs[seg_idx].dir ^ invert;
1472               art_svp_add_segment (&new_vp, &n_segs_max,
1473                                    NULL,
1474                                    n_points, new_dir, new_points,
1475                                    &vp->segs[seg_idx].bbox);
1476             }
1477
1478           tmp1 = seg_idx;
1479           for (j = i; j < n_active_segs; j++)
1480             {
1481               tmp2 = active_segs[j];
1482               active_segs[j] = tmp1;
1483               tmp1 = tmp2;
1484             }
1485           active_segs[n_active_segs] = tmp1;
1486           n_active_segs++;
1487           seg_idx++;
1488         }
1489
1490 #ifdef VERBOSE
1491       /* all active segs cross the y scanline (considering segs to be
1492        closed on top and open on bottom) */
1493       for (i = 0; i < n_active_segs; i++)
1494         {
1495           asi = active_segs[i];
1496           printf ("%d:%d (%g, %g) - (%g, %g) %s %d\n", asi,
1497                   cursor[asi],
1498                   vp->segs[asi].points[cursor[asi]].x,
1499                   vp->segs[asi].points[cursor[asi]].y,
1500                   vp->segs[asi].points[cursor[asi] + 1].x,
1501                   vp->segs[asi].points[cursor[asi] + 1].y,
1502                   vp->segs[asi].dir ? "v" : "^",
1503                   winding[asi]);
1504         }
1505 #endif
1506
1507       /* advance y to the next event */
1508       if (n_active_segs == 0)
1509         {
1510           if (seg_idx < vp->n_segs)
1511             y = vp->segs[seg_idx].points[0].y;
1512           /* else we're done */
1513         }
1514       else
1515         {
1516           asi = active_segs[0];
1517           y = vp->segs[asi].points[cursor[asi] + 1].y;
1518           for (i = 1; i < n_active_segs; i++)
1519             {
1520               asi = active_segs[i];
1521               if (y > vp->segs[asi].points[cursor[asi] + 1].y)
1522                 y = vp->segs[asi].points[cursor[asi] + 1].y;
1523             }
1524           if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
1525             y = vp->segs[seg_idx].points[0].y;
1526         }
1527
1528       /* advance cursors to reach new y */
1529       for (i = 0; i < n_active_segs; i++)
1530         {
1531           asi = active_segs[i];
1532           while (cursor[asi] < vp->segs[asi].n_points - 1 &&
1533                  y >= vp->segs[asi].points[cursor[asi] + 1].y)
1534             cursor[asi]++;
1535         }
1536 #ifdef VERBOSE
1537       printf ("\n");
1538 #endif
1539     }
1540   art_free (cursor);
1541   art_free (active_segs);
1542   art_free (winding);
1543
1544   return new_vp;
1545 }