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