1 /* Libart_LGPL - library of basic graphic primitives
2 * Copyright (C) 1998-2000 Raph Levien
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.
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.
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.
20 /* Primitive intersection and winding number operations on sorted
23 These routines are internal to libart, used to construct operations
24 like intersection, union, and difference. */
27 #include <stdio.h> /* for printf of debugging info */
28 #include <string.h> /* for memcpy */
34 #include "art_svp_wind.h"
38 #define PT_EQ(p1,p2) ((p1).x == (p2).x && (p1).y == (p2).y)
40 #define PT_CLOSE(p1,p2) (fabs ((p1).x - (p2).x) < 1e-6 && fabs ((p1).y - (p2).y) < 1e-6)
42 /* return nonzero and set *p to the intersection point if the lines
43 z0-z1 and z2-z3 intersect each other. */
45 intersect_lines (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3,
50 double d0, d1, d2, d3;
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))
58 if (PT_CLOSE (z0, z2) || PT_CLOSE (z0, z3) || PT_CLOSE (z1, z2) || PT_CLOSE (z1, z3))
62 /* find line equations ax + by + c = 0 */
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) */
69 d2 = a01 * z2.x + b01 * z2.y + c01;
70 d3 = a01 * z3.x + b01 * z3.y + c01;
71 if ((d2 > 0) == (d3 > 0))
76 c23 = -(z2.x * a23 + z2.y * b23);
78 d0 = a23 * z0.x + b23 * z0.y + c23;
79 d1 = a23 * z1.x + b23 * z1.y + c23;
80 if ((d0 > 0) == (d1 > 0))
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);
95 trap_epsilon (double v)
97 const double epsilon = EPSILON;
99 if (v < epsilon && v > -epsilon) return 0;
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,
108 The case analysis in this function is quite ugly. The fact that it's
109 almost 200 lines long is ridiculous.
111 Ok, so here's the plan to cut it down:
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.
117 Then, do d0 and d1 computation, but only if a23 is nonzero.
119 Finally, do d2 and d3 computation, but only if a01 is nonzero.
121 Fall through to returning 0 (this will happen when both lines are
122 horizontal and they overlap).
125 x_order (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3)
127 double a01, b01, c01;
128 double a23, b23, c23;
129 double d0, d1, d2, d3;
135 double x01min, x01max;
136 double x23min, x23max;
160 if (x23min >= x01max) return 1;
161 else if (x01min >= x23max) return -1;
166 /* z0-z1 is horizontal, z2-z3 isn't */
169 c23 = -(z2.x * a23 + z2.y * b23);
178 d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23);
179 d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23);
183 if (d1 >= 0) return 1;
188 if (d1 > 0) return 1;
189 else if (d1 < 0) return -1;
190 else printf ("case 1 degenerate\n");
195 if (d1 <= 0) return -1;
200 else if (z2.y == z3.y)
202 /* z2-z3 is horizontal, z0-z1 isn't */
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) */
216 d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01);
217 d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01);
221 if (d3 >= 0) return -1;
226 if (d3 > 0) return -1;
227 else if (d3 < 0) return 1;
228 else printf ("case 2 degenerate\n");
233 if (d3 <= 0) return 1;
238 /* find line equations ax + by + c = 0 */
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) */
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 */
254 d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01);
255 d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01);
258 if (d3 >= 0) return -1;
262 if (d3 > 0) return -1;
263 else if (d3 < 0) return 1;
265 fprintf (stderr, "colinear!\n");
269 if (d3 <= 0) return 1;
274 c23 = -(z2.x * a23 + z2.y * b23);
282 d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23);
283 d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23);
286 if (d1 >= 0) return 1;
290 if (d1 > 0) return 1;
291 else if (d1 < 0) return -1;
293 fprintf (stderr, "colinear!\n");
297 if (d1 <= 0) return -1;
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 */
306 x_order_2 (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3)
308 double a23, b23, c23;
313 c23 = -(z2.x * a23 + z2.y * b23);
322 d0 = a23 * z0.x + b23 * z0.y + c23;
326 else if (d0 < -EPSILON)
329 d1 = a23 * z1.x + b23 * z1.y + c23;
332 else if (d1 < -EPSILON)
335 if (z0.x <= z2.x && z1.x <= z2.x && z0.x <= z3.x && z1.x <= z3.x)
337 if (z0.x >= z2.x && z1.x >= z2.x && z0.x >= z3.x && z1.x >= z3.x)
340 fprintf (stderr, "x_order_2: colinear!\n");
345 /* Traverse the vector path, keeping it in x-sorted order.
347 This routine doesn't actually do anything - it's just here for
348 explanatory purposes. */
350 traverse (ArtSVP *vp)
361 active_segs = art_new (int, vp->n_segs);
362 cursor = art_new (int, vp->n_segs);
366 y = vp->segs[0].points[0].y;
367 while (seg_idx < vp->n_segs || n_active_segs > 0)
369 printf ("y = %g\n", y);
370 /* delete segments ending at y from active list */
371 for (i = 0; i < n_active_segs; i++)
373 asi = active_segs[i];
374 if (vp->segs[asi].n_points - 1 == cursor[asi] &&
375 vp->segs[asi].points[cursor[asi]].y == y)
377 printf ("deleting %d\n", asi);
379 for (j = i; j < n_active_segs; j++)
380 active_segs[j] = active_segs[j + 1];
385 /* insert new segments into the active list */
386 while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
389 printf ("inserting %d\n", seg_idx);
390 for (i = 0; i < n_active_segs; i++)
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)
400 for (j = i; j < n_active_segs; j++)
402 tmp2 = active_segs[j];
403 active_segs[j] = tmp1;
406 active_segs[n_active_segs] = tmp1;
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++)
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" : "^");
424 /* advance y to the next event */
425 if (n_active_segs == 0)
427 if (seg_idx < vp->n_segs)
428 y = vp->segs[seg_idx].points[0].y;
429 /* else we're done */
433 asi = active_segs[0];
434 y = vp->segs[asi].points[cursor[asi] + 1].y;
435 for (i = 1; i < n_active_segs; i++)
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;
441 if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
442 y = vp->segs[seg_idx].points[0].y;
445 /* advance cursors to reach new y */
446 for (i = 0; i < n_active_segs; i++)
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)
456 art_free (active_segs);
460 /* I believe that the loop will always break with i=1.
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
469 This change should cut down on allocation bandwidth, and also
470 eliminate the iteration through n_ipl below.
474 insert_ip (int seg_i, int *n_ips, int *n_ips_max, ArtPoint **ips, ArtPoint ip)
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]);
485 for (i = 1; i < n_ipl; i++)
489 for (; i <= n_ipl; i++)
497 /* test active segment (i - 1) against i for intersection, if
498 so, add intersection point to both ips lists. */
500 intersect_neighbors (int i, int *active_segs,
501 int *n_ips, int *n_ips_max, ArtPoint **ips,
502 int *cursor, ArtSVP *vp)
504 ArtPoint z0, z1, z2, z3;
508 asi01 = active_segs[i - 1];
511 if (n_ips[asi01] == 1)
512 z1 = vp->segs[asi01].points[cursor[asi01] + 1];
516 asi23 = active_segs[i];
519 if (n_ips[asi23] == 1)
520 z3 = vp->segs[asi23].points[cursor[asi23] + 1];
524 if (intersect_lines (z0, z1, z2, z3, &ip))
527 printf ("new intersection point: (%g, %g)\n", ip.x, ip.y);
529 insert_ip (asi01, n_ips, n_ips_max, ips, ip);
530 insert_ip (asi23, n_ips, n_ips_max, ips, ip);
534 /* Add a new point to a segment in the svp.
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
541 svp_add_point (ArtSVP *svp, int *n_points_max,
542 ArtPoint p, int *seg_map, int *active_segs, int n_active_segs,
545 int asi, asi_left, asi_right;
546 int n_points, n_points_left, n_points_right;
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 */
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]))
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],
570 printf ("svp_add_point: cross on left!\n");
576 if (i + 1 < n_active_segs)
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]))
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],
592 printf ("svp_add_point: cross on right!\n");
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)
602 else if (p.x > seg->bbox.x1)
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.
612 "Where it is supposed to go" is defined as following:
614 Delete element i, re-insert at position target (bumping everything
615 target and greater to the right).
618 find_crossing (int i, int *active_segs, int n_active_segs,
619 int *cursor, ArtPoint **ips, int *n_ips, ArtSVP *vp)
621 int asi, asi_left, asi_right;
627 asi = active_segs[i];
630 p1 = vp->segs[asi].points[cursor[asi] + 1];
634 for (target = i; target > 0; target--)
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];
641 p1l = ips[asi_left][1];
642 if (!PT_EQ (p0, p0l))
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);
649 if (x_order (p0l, p1l, p0, p1) == 1)
653 printf ("scanning to the left (i=%d, target=%d)\n", i, target);
657 if (target < i) return target;
659 for (; target < n_active_segs - 1; target++)
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];
666 p1r = ips[asi_right][1];
667 if (!PT_EQ (p0, p0r))
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);
674 if (x_order (p0r, p1r, p0, p1) == 1)
678 printf ("scanning to the right (i=%d, target=%d)\n", i, target);
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.
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.
694 This routine takes care of intersecting the appropriate neighbors,
697 A first argument of -1 immediately returns, which helps reduce special
698 casing in the main unwind routine.
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,
719 printf ("fix_crossing: [%d..%d)", start, end);
720 for (k = 0; k < n_active_segs; k++)
721 printf (" %d", active_segs[k]);
728 for (i = start + 1; i < end; i++)
731 asi = active_segs[i];
732 if (cursor[asi] < vp->segs[asi].n_points - 1) {
735 p1i = vp->segs[asi].points[cursor[asi] + 1];
739 for (j = i - 1; j >= start; j--)
741 asj = active_segs[j];
742 if (cursor[asj] < vp->segs[asj].n_points - 1)
746 p1j = vp->segs[asj].points[cursor[asj] + 1];
750 /* we _hope_ p0i = p0j */
751 if (x_order_2 (p0j, p1j, p0i, p1i) == -1)
757 /* target is where active_seg[i] _should_ be in active_segs */
764 printf ("fix_crossing: at %i should be %i\n", i, target);
767 /* let's close off all relevant segments */
768 for (j = i; j >= target; j--)
770 asi = active_segs[j];
771 /* First conjunct: this isn't the last point in the original
774 Second conjunct: this isn't the first point in the new
775 segment (i.e. already broken).
777 if (cursor[asi] < vp->segs[asi].n_points - 1 &&
778 (*p_new_vp)->segs[seg_map[asi]].n_points != 1)
783 printf ("closing off %d\n", j);
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,
790 1, vp->segs[asi].dir,
793 (*pn_points_max)[seg_num] = 16;
794 seg_map[asi] = seg_num;
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;
806 if (swap && start > 0)
810 as_start = active_segs[start];
811 if (cursor[as_start] < vp->segs[as_start].n_points)
814 printf ("checking intersection of %d, %d\n", start - 1, start);
816 intersect_neighbors (start, active_segs,
817 n_ips, n_ips_max, ips,
822 if (swap && end < n_active_segs)
826 as_end = active_segs[end - 1];
827 if (cursor[as_end] < vp->segs[as_end].n_points)
830 printf ("checking intersection of %d, %d\n", end - 1, end);
832 intersect_neighbors (end, active_segs,
833 n_ips, n_ips_max, ips,
840 printf ("fix_crossing return: [%d..%d)", start, end);
841 for (k = 0; k < n_active_segs; k++)
842 printf (" %d", active_segs[k]);
848 /* Return a new sorted vector that covers the same area as the
849 argument, but which satisfies the nocross invariant.
851 Basically, this routine works by finding the intersection points,
852 and cutting the segments at those points.
854 Status of this routine:
856 Basic correctness: Seems ok.
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.
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.
866 Precision: pretty good, although the numerical stability problems
867 make this routine unsuitable for precise calculations of
872 /* Here is a more detailed description of the algorithm. It follows
873 roughly the structure of traverse (above), but is obviously quite
876 Here are a few important data structures:
878 A new sorted vector path (new_svp).
880 For each (active) segment in the original, a list of intersection
883 Of course, the original being traversed.
885 The following invariants hold (in addition to the invariants
886 of the traverse procedure).
888 The new sorted vector path lies entirely above the y scan line.
890 The new sorted vector path keeps the nocross invariant.
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).
896 The list of intersection points + the (cursor + 1) point is kept
897 in nondecreasing y order.
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).
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
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.
915 Another important data structure is mapping from original segment
916 number to new segment number.
918 The algorithm is perhaps best understood as advancing the cursors
919 while maintaining these invariants. Here's roughly how it's done.
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).
925 When inserting new segments, they are intersection tested against
926 their neighbors. The top point of the segment becomes the first
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.
937 Two steps need further clarification:
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).
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.
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. */
959 * art_svp_uncross: Resolve self-intersections of an svp.
960 * @vp: The original svp.
962 * Finds all the intersections within @vp, and constructs a new svp
963 * with new points added at these intersections.
965 * This routine needs to be redone from scratch with numerical robustness
966 * in mind. I'm working on it.
968 * Return value: The new svp.
971 art_svp_uncross (ArtSVP *vp)
981 /* new data structures */
982 /* intersection points; invariant: *ips[i] is only allocated if
984 int *n_ips, *n_ips_max;
986 /* new sorted vector path */
987 int n_segs_max, seg_num;
990 /* mapping from argument to new segment numbers - again, only valid
1000 new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) +
1001 (n_segs_max - 1) * sizeof(ArtSVPSeg));
1004 if (vp->n_segs == 0)
1007 active_segs = art_new (int, vp->n_segs);
1008 cursor = art_new (int, vp->n_segs);
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);
1015 n_points_max = art_new (int, n_segs_max);
1019 y = vp->segs[0].points[0].y;
1020 while (seg_idx < vp->n_segs || n_active_segs > 0)
1023 printf ("y = %g\n", y);
1026 /* maybe move deletions to end of loop (to avoid so much special
1027 casing on the end of a segment)? */
1029 /* delete segments ending at y from active list */
1030 for (i = 0; i < n_active_segs; i++)
1032 asi = active_segs[i];
1033 if (vp->segs[asi].n_points - 1 == cursor[asi] &&
1034 vp->segs[asi].points[cursor[asi]].y == y)
1039 printf ("deleting %d\n", asi);
1041 art_free (ips[asi]);
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];
1050 while (vp->segs[asi].n_points - 1 == cursor[asi] &&
1051 vp->segs[asi].points[cursor[asi]].y == y);
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,
1063 /* insert new segments into the active list */
1064 while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
1067 printf ("inserting %d\n", seg_idx);
1069 cursor[seg_idx] = 0;
1070 for (i = 0; i < n_active_segs; i++)
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)
1080 /* Create and initialize the intersection points data structure */
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];
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,
1091 1, vp->segs[seg_idx].dir,
1094 n_points_max[seg_num] = 16;
1095 seg_map[seg_idx] = seg_num;
1098 for (j = i; j < n_active_segs; j++)
1100 tmp2 = active_segs[j];
1101 active_segs[j] = tmp1;
1104 active_segs[n_active_segs] = tmp1;
1108 intersect_neighbors (i, active_segs,
1109 n_ips, n_ips_max, ips,
1112 if (i + 1 < n_active_segs)
1113 intersect_neighbors (i + 1, active_segs,
1114 n_ips, n_ips_max, ips,
1120 /* all active segs cross the y scanline (considering segs to be
1121 closed on top and open on bottom) */
1123 for (i = 0; i < n_active_segs; i++)
1125 asi = active_segs[i];
1126 printf ("%d ", asi);
1127 for (j = 0; j < n_ips[asi]; j++)
1128 printf ("(%g, %g) - ",
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" : "^");
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)
1143 if (seg_idx < vp->n_segs)
1144 y = vp->segs[seg_idx].points[0].y;
1145 /* else we're done */
1149 asi = active_segs[0];
1150 if (n_ips[asi] == 1)
1151 y = vp->segs[asi].points[cursor[asi] + 1].y;
1154 for (i = 1; i < n_active_segs; i++)
1156 asi = active_segs[i];
1157 if (n_ips[asi] == 1)
1158 y_curs = vp->segs[asi].points[cursor[asi] + 1].y;
1160 y_curs = ips[asi][1].y;
1164 if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
1165 y = vp->segs[seg_idx].points[0].y;
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++)
1174 asi = active_segs[i];
1175 if (n_ips[asi] == 1)
1176 p_curs = vp->segs[asi].points[cursor[asi] + 1];
1178 p_curs = ips[asi][1];
1181 svp_add_point (new_vp, n_points_max,
1182 p_curs, seg_map, active_segs, n_active_segs, i);
1185 for (j = 0; j < n_ips[asi]; j++)
1186 ips[asi][j] = ips[asi][j + 1];
1188 if (n_ips[asi] == 0)
1190 ips[asi][0] = p_curs;
1195 if (first_share < 0 || p_curs.x != share_x)
1197 /* this is where crossings are detected, and if
1198 found, the active segments switched around. */
1200 fix_crossing (first_share, i,
1201 active_segs, n_active_segs,
1202 cursor, ips, n_ips, n_ips_max, vp, seg_map,
1204 &n_segs_max, &n_points_max);
1210 if (cursor[asi] < vp->segs[asi].n_points - 1)
1214 intersect_neighbors (i, active_segs,
1215 n_ips, n_ips_max, ips,
1218 if (i + 1 < n_active_segs)
1219 intersect_neighbors (i + 1, active_segs,
1220 n_ips, n_ips_max, ips,
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,
1231 &n_segs_max, &n_points_max);
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,
1241 &n_segs_max, &n_points_max);
1248 /* not necessary to sort, new segments only get added at y, which
1249 increases monotonically */
1251 qsort (&new_vp->segs, new_vp->n_segs, sizeof (svp_seg), svp_seg_compare);
1254 for (k = 0; k < new_vp->n_segs - 1; k++)
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);
1270 art_free (n_points_max);
1272 art_free (n_ips_max);
1276 art_free (active_segs);
1283 /* Rewind a svp satisfying the nocross invariant.
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.
1290 Status of this routine:
1292 Basic correctness: Was ok in gfonted. However, this code does not
1293 yet compute bboxes for the resulting svp segs.
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.
1301 Precision: good, except that no attempt is made to remove "hair".
1302 Doing random perturbation just makes matters worse.
1306 * art_svp_rewind_uncrossed: Rewind an svp satisfying the nocross invariant.
1307 * @vp: The original svp.
1308 * @rule: The winding rule.
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.
1314 * With @rule == ART_WIND_RULE_NONZERO, the resulting svp has a
1315 * winding number of 1 where @vp has a nonzero winding number.
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.
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.
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.
1330 * This routine needs to be redone from scratch with numerical robustness
1331 * in mind. I'm working on it.
1333 * Return value: The new svp.
1336 art_svp_rewind_uncrossed (ArtSVP *vp, ArtWindRule rule)
1358 new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) +
1359 (n_segs_max - 1) * sizeof(ArtSVPSeg));
1362 if (vp->n_segs == 0)
1365 winding = art_new (int, vp->n_segs);
1367 active_segs = art_new (int, vp->n_segs);
1368 cursor = art_new (int, vp->n_segs);
1372 y = vp->segs[0].points[0].y;
1373 while (seg_idx < vp->n_segs || n_active_segs > 0)
1376 printf ("y = %g\n", y);
1378 /* delete segments ending at y from active list */
1379 for (i = 0; i < n_active_segs; i++)
1381 asi = active_segs[i];
1382 if (vp->segs[asi].n_points - 1 == cursor[asi] &&
1383 vp->segs[asi].points[cursor[asi]].y == y)
1386 printf ("deleting %d\n", asi);
1389 for (j = i; j < n_active_segs; j++)
1390 active_segs[j] = active_segs[j + 1];
1395 /* insert new segments into the active list */
1396 while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
1399 printf ("inserting %d\n", seg_idx);
1401 cursor[seg_idx] = 0;
1402 for (i = 0; i < n_active_segs; i++)
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)
1412 /* Determine winding number for this segment */
1415 else if (vp->segs[active_segs[i - 1]].dir)
1416 left_wind = winding[active_segs[i - 1]];
1418 left_wind = winding[active_segs[i - 1]] - 1;
1420 if (vp->segs[seg_idx].dir)
1421 wind = left_wind + 1;
1425 winding[seg_idx] = wind;
1429 case ART_WIND_RULE_NONZERO:
1430 keep = (wind == 1 || wind == 0);
1431 invert = (wind == 0);
1433 case ART_WIND_RULE_INTERSECT:
1437 case ART_WIND_RULE_ODDEVEN:
1439 invert = !(wind & 1);
1441 case ART_WIND_RULE_POSITIVE:
1453 ArtPoint *points, *new_points;
1458 printf ("keeping segment %d\n", seg_idx);
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,
1467 n_points, new_dir, new_points,
1468 &vp->segs[seg_idx].bbox);
1472 for (j = i; j < n_active_segs; j++)
1474 tmp2 = active_segs[j];
1475 active_segs[j] = tmp1;
1478 active_segs[n_active_segs] = tmp1;
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++)
1488 asi = active_segs[i];
1489 printf ("%d:%d (%g, %g) - (%g, %g) %s %d\n", 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" : "^",
1500 /* advance y to the next event */
1501 if (n_active_segs == 0)
1503 if (seg_idx < vp->n_segs)
1504 y = vp->segs[seg_idx].points[0].y;
1505 /* else we're done */
1509 asi = active_segs[0];
1510 y = vp->segs[asi].points[cursor[asi] + 1].y;
1511 for (i = 1; i < n_active_segs; i++)
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;
1517 if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
1518 y = vp->segs[seg_idx].points[0].y;
1521 /* advance cursors to reach new y */
1522 for (i = 0; i < n_active_segs; i++)
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)
1534 art_free (active_segs);