added implementation for colortransform
[swftools.git] / 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 }