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 "art_svp_wind.h"
29 #include <stdio.h> /* for printf of debugging info */
30 #include <string.h> /* for memcpy */
39 #define PT_EQ(p1,p2) ((p1).x == (p2).x && (p1).y == (p2).y)
41 #define PT_CLOSE(p1,p2) (fabs ((p1).x - (p2).x) < 1e-6 && fabs ((p1).y - (p2).y) < 1e-6)
43 /* return nonzero and set *p to the intersection point if the lines
44 z0-z1 and z2-z3 intersect each other. */
46 intersect_lines (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3,
51 double d0, d1, d2, d3;
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))
59 if (PT_CLOSE (z0, z2) || PT_CLOSE (z0, z3) || PT_CLOSE (z1, z2) || PT_CLOSE (z1, z3))
63 /* find line equations ax + by + c = 0 */
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) */
70 d2 = a01 * z2.x + b01 * z2.y + c01;
71 d3 = a01 * z3.x + b01 * z3.y + c01;
72 if ((d2 > 0) == (d3 > 0))
77 c23 = -(z2.x * a23 + z2.y * b23);
79 d0 = a23 * z0.x + b23 * z0.y + c23;
80 d1 = a23 * z1.x + b23 * z1.y + c23;
81 if ((d0 > 0) == (d1 > 0))
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);
96 trap_epsilon (double v)
98 const double epsilon = EPSILON;
100 if (v < epsilon && v > -epsilon) return 0;
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,
109 The case analysis in this function is quite ugly. The fact that it's
110 almost 200 lines long is ridiculous.
112 Ok, so here's the plan to cut it down:
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.
118 Then, do d0 and d1 computation, but only if a23 is nonzero.
120 Finally, do d2 and d3 computation, but only if a01 is nonzero.
122 Fall through to returning 0 (this will happen when both lines are
123 horizontal and they overlap).
126 x_order (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3)
128 double a01, b01, c01;
129 double a23, b23, c23;
130 double d0, d1, d2, d3;
136 double x01min, x01max;
137 double x23min, x23max;
161 if (x23min >= x01max) return 1;
162 else if (x01min >= x23max) return -1;
167 /* z0-z1 is horizontal, z2-z3 isn't */
170 c23 = -(z2.x * a23 + z2.y * b23);
179 d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23);
180 d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23);
184 if (d1 >= 0) return 1;
189 if (d1 > 0) return 1;
190 else if (d1 < 0) return -1;
191 else /*printf ("case 1 degenerate\n")*/;
196 if (d1 <= 0) return -1;
201 else if (z2.y == z3.y)
203 /* z2-z3 is horizontal, z0-z1 isn't */
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) */
217 d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01);
218 d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01);
222 if (d3 >= 0) return -1;
227 if (d3 > 0) return -1;
228 else if (d3 < 0) return 1;
229 else printf ("case 2 degenerate\n");
234 if (d3 <= 0) return 1;
239 /* find line equations ax + by + c = 0 */
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) */
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 */
255 d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01);
256 d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01);
259 if (d3 >= 0) return -1;
263 if (d3 > 0) return -1;
264 else if (d3 < 0) return 1;
266 /*fprintf (stderr, "colinear!\n")*/;
270 if (d3 <= 0) return 1;
275 c23 = -(z2.x * a23 + z2.y * b23);
283 d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23);
284 d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23);
287 if (d1 >= 0) return 1;
291 if (d1 > 0) return 1;
292 else if (d1 < 0) return -1;
294 /*fprintf (stderr, "colinear!\n")*/;
298 if (d1 <= 0) return -1;
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 */
307 x_order_2 (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3)
309 double a23, b23, c23;
314 c23 = -(z2.x * a23 + z2.y * b23);
323 d0 = a23 * z0.x + b23 * z0.y + c23;
327 else if (d0 < -EPSILON)
330 d1 = a23 * z1.x + b23 * z1.y + c23;
333 else if (d1 < -EPSILON)
336 if (z0.x == z1.x && z1.x == z2.x && z2.x == z3.x)
338 //art_dprint ("x_order_2: colinear and horizontally aligned!\n");
342 if (z0.x <= z2.x && z1.x <= z2.x && z0.x <= z3.x && z1.x <= z3.x)
344 if (z0.x >= z2.x && z1.x >= z2.x && z0.x >= z3.x && z1.x >= z3.x)
347 //fprintf (stderr, "x_order_2: colinear!\n");
352 /* Traverse the vector path, keeping it in x-sorted order.
354 This routine doesn't actually do anything - it's just here for
355 explanatory purposes. */
357 traverse (ArtSVP *vp)
368 active_segs = art_new (int, vp->n_segs);
369 cursor = art_new (int, vp->n_segs);
373 y = vp->segs[0].points[0].y;
374 while (seg_idx < vp->n_segs || n_active_segs > 0)
376 printf ("y = %g\n", y);
377 /* delete segments ending at y from active list */
378 for (i = 0; i < n_active_segs; i++)
380 asi = active_segs[i];
381 if (vp->segs[asi].n_points - 1 == cursor[asi] &&
382 vp->segs[asi].points[cursor[asi]].y == y)
384 printf ("deleting %d\n", asi);
386 for (j = i; j < n_active_segs; j++)
387 active_segs[j] = active_segs[j + 1];
392 /* insert new segments into the active list */
393 while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
396 printf ("inserting %d\n", seg_idx);
397 for (i = 0; i < n_active_segs; i++)
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)
407 for (j = i; j < n_active_segs; j++)
409 tmp2 = active_segs[j];
410 active_segs[j] = tmp1;
413 active_segs[n_active_segs] = tmp1;
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++)
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" : "^");
431 /* advance y to the next event */
432 if (n_active_segs == 0)
434 if (seg_idx < vp->n_segs)
435 y = vp->segs[seg_idx].points[0].y;
436 /* else we're done */
440 asi = active_segs[0];
441 y = vp->segs[asi].points[cursor[asi] + 1].y;
442 for (i = 1; i < n_active_segs; i++)
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;
448 if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
449 y = vp->segs[seg_idx].points[0].y;
452 /* advance cursors to reach new y */
453 for (i = 0; i < n_active_segs; i++)
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)
463 art_free (active_segs);
467 /* I believe that the loop will always break with i=1.
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
476 This change should cut down on allocation bandwidth, and also
477 eliminate the iteration through n_ipl below.
481 insert_ip (int seg_i, int *n_ips, int *n_ips_max, ArtPoint **ips, ArtPoint ip)
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]);
492 for (i = 1; i < n_ipl; i++)
496 for (; i <= n_ipl; i++)
504 /* test active segment (i - 1) against i for intersection, if
505 so, add intersection point to both ips lists. */
507 intersect_neighbors (int i, int *active_segs,
508 int *n_ips, int *n_ips_max, ArtPoint **ips,
509 int *cursor, ArtSVP *vp)
511 ArtPoint z0, z1, z2, z3;
515 asi01 = active_segs[i - 1];
518 if (n_ips[asi01] == 1)
519 z1 = vp->segs[asi01].points[cursor[asi01] + 1];
523 asi23 = active_segs[i];
526 if (n_ips[asi23] == 1)
527 z3 = vp->segs[asi23].points[cursor[asi23] + 1];
531 if (intersect_lines (z0, z1, z2, z3, &ip))
534 printf ("new intersection point: (%g, %g)\n", ip.x, ip.y);
536 insert_ip (asi01, n_ips, n_ips_max, ips, ip);
537 insert_ip (asi23, n_ips, n_ips_max, ips, ip);
541 /* Add a new point to a segment in the svp.
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
548 svp_add_point (ArtSVP *svp, int *n_points_max,
549 ArtPoint p, int *seg_map, int *active_segs, int n_active_segs,
552 int asi, asi_left, asi_right;
553 int n_points, n_points_left, n_points_right;
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 */
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]))
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],
577 printf ("svp_add_point: cross on left!\n");
583 if (i + 1 < n_active_segs)
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]))
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],
599 printf ("svp_add_point: cross on right!\n");
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)
609 else if (p.x > seg->bbox.x1)
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.
619 "Where it is supposed to go" is defined as following:
621 Delete element i, re-insert at position target (bumping everything
622 target and greater to the right).
625 find_crossing (int i, int *active_segs, int n_active_segs,
626 int *cursor, ArtPoint **ips, int *n_ips, ArtSVP *vp)
628 int asi, asi_left, asi_right;
634 asi = active_segs[i];
637 p1 = vp->segs[asi].points[cursor[asi] + 1];
641 for (target = i; target > 0; target--)
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];
648 p1l = ips[asi_left][1];
649 if (!PT_EQ (p0, p0l))
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);
656 if (x_order (p0l, p1l, p0, p1) == 1)
660 printf ("scanning to the left (i=%d, target=%d)\n", i, target);
664 if (target < i) return target;
666 for (; target < n_active_segs - 1; target++)
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];
673 p1r = ips[asi_right][1];
674 if (!PT_EQ (p0, p0r))
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);
681 if (x_order (p0r, p1r, p0, p1) == 1)
685 printf ("scanning to the right (i=%d, target=%d)\n", i, target);
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.
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.
701 This routine takes care of intersecting the appropriate neighbors,
704 A first argument of -1 immediately returns, which helps reduce special
705 casing in the main unwind routine.
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,
726 printf ("fix_crossing: [%d..%d)", start, end);
727 for (k = 0; k < n_active_segs; k++)
728 printf (" %d", active_segs[k]);
735 for (i = start + 1; i < end; i++)
738 asi = active_segs[i];
739 if (cursor[asi] < vp->segs[asi].n_points - 1) {
742 p1i = vp->segs[asi].points[cursor[asi] + 1];
746 for (j = i - 1; j >= start; j--)
748 asj = active_segs[j];
749 if (cursor[asj] < vp->segs[asj].n_points - 1)
753 p1j = vp->segs[asj].points[cursor[asj] + 1];
757 /* we _hope_ p0i = p0j */
758 if (x_order_2 (p0j, p1j, p0i, p1i) == -1)
764 /* target is where active_seg[i] _should_ be in active_segs */
771 printf ("fix_crossing: at %i should be %i\n", i, target);
774 /* let's close off all relevant segments */
775 for (j = i; j >= target; j--)
777 asi = active_segs[j];
778 /* First conjunct: this isn't the last point in the original
781 Second conjunct: this isn't the first point in the new
782 segment (i.e. already broken).
784 if (cursor[asi] < vp->segs[asi].n_points - 1 &&
785 (*p_new_vp)->segs[seg_map[asi]].n_points != 1)
790 printf ("closing off %d\n", j);
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,
797 1, vp->segs[asi].dir,
800 (*pn_points_max)[seg_num] = 16;
801 seg_map[asi] = seg_num;
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;
813 if (swap && start > 0)
817 as_start = active_segs[start];
818 if (cursor[as_start] < vp->segs[as_start].n_points)
821 printf ("checking intersection of %d, %d\n", start - 1, start);
823 intersect_neighbors (start, active_segs,
824 n_ips, n_ips_max, ips,
829 if (swap && end < n_active_segs)
833 as_end = active_segs[end - 1];
834 if (cursor[as_end] < vp->segs[as_end].n_points)
837 printf ("checking intersection of %d, %d\n", end - 1, end);
839 intersect_neighbors (end, active_segs,
840 n_ips, n_ips_max, ips,
847 printf ("fix_crossing return: [%d..%d)", start, end);
848 for (k = 0; k < n_active_segs; k++)
849 printf (" %d", active_segs[k]);
855 /* Return a new sorted vector that covers the same area as the
856 argument, but which satisfies the nocross invariant.
858 Basically, this routine works by finding the intersection points,
859 and cutting the segments at those points.
861 Status of this routine:
863 Basic correctness: Seems ok.
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.
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.
873 Precision: pretty good, although the numerical stability problems
874 make this routine unsuitable for precise calculations of
879 /* Here is a more detailed description of the algorithm. It follows
880 roughly the structure of traverse (above), but is obviously quite
883 Here are a few important data structures:
885 A new sorted vector path (new_svp).
887 For each (active) segment in the original, a list of intersection
890 Of course, the original being traversed.
892 The following invariants hold (in addition to the invariants
893 of the traverse procedure).
895 The new sorted vector path lies entirely above the y scan line.
897 The new sorted vector path keeps the nocross invariant.
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).
903 The list of intersection points + the (cursor + 1) point is kept
904 in nondecreasing y order.
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).
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
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.
922 Another important data structure is mapping from original segment
923 number to new segment number.
925 The algorithm is perhaps best understood as advancing the cursors
926 while maintaining these invariants. Here's roughly how it's done.
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).
932 When inserting new segments, they are intersection tested against
933 their neighbors. The top point of the segment becomes the first
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.
944 Two steps need further clarification:
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).
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.
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. */
966 * art_svp_uncross: Resolve self-intersections of an svp.
967 * @vp: The original svp.
969 * Finds all the intersections within @vp, and constructs a new svp
970 * with new points added at these intersections.
972 * This routine needs to be redone from scratch with numerical robustness
973 * in mind. I'm working on it.
975 * Return value: The new svp.
978 art_svp_uncross (ArtSVP *vp)
988 /* new data structures */
989 /* intersection points; invariant: *ips[i] is only allocated if
991 int *n_ips, *n_ips_max;
993 /* new sorted vector path */
994 int n_segs_max, seg_num;
997 /* mapping from argument to new segment numbers - again, only valid
1007 new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) +
1008 (n_segs_max - 1) * sizeof(ArtSVPSeg));
1011 if (vp->n_segs == 0)
1014 active_segs = art_new (int, vp->n_segs);
1015 cursor = art_new (int, vp->n_segs);
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);
1022 n_points_max = art_new (int, n_segs_max);
1026 y = vp->segs[0].points[0].y;
1027 while (seg_idx < vp->n_segs || n_active_segs > 0)
1030 printf ("y = %g\n", y);
1033 /* maybe move deletions to end of loop (to avoid so much special
1034 casing on the end of a segment)? */
1036 /* delete segments ending at y from active list */
1037 for (i = 0; i < n_active_segs; i++)
1039 asi = active_segs[i];
1040 if (vp->segs[asi].n_points - 1 == cursor[asi] &&
1041 vp->segs[asi].points[cursor[asi]].y == y)
1046 printf ("deleting %d\n", asi);
1048 art_free (ips[asi]);
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];
1057 while (vp->segs[asi].n_points - 1 == cursor[asi] &&
1058 vp->segs[asi].points[cursor[asi]].y == y);
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,
1070 /* insert new segments into the active list */
1071 while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
1074 printf ("inserting %d\n", seg_idx);
1076 cursor[seg_idx] = 0;
1077 for (i = 0; i < n_active_segs; i++)
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)
1087 /* Create and initialize the intersection points data structure */
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];
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,
1098 1, vp->segs[seg_idx].dir,
1101 n_points_max[seg_num] = 16;
1102 seg_map[seg_idx] = seg_num;
1105 for (j = i; j < n_active_segs; j++)
1107 tmp2 = active_segs[j];
1108 active_segs[j] = tmp1;
1111 active_segs[n_active_segs] = tmp1;
1115 intersect_neighbors (i, active_segs,
1116 n_ips, n_ips_max, ips,
1119 if (i + 1 < n_active_segs)
1120 intersect_neighbors (i + 1, active_segs,
1121 n_ips, n_ips_max, ips,
1127 /* all active segs cross the y scanline (considering segs to be
1128 closed on top and open on bottom) */
1130 for (i = 0; i < n_active_segs; i++)
1132 asi = active_segs[i];
1133 printf ("%d ", asi);
1134 for (j = 0; j < n_ips[asi]; j++)
1135 printf ("(%g, %g) - ",
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" : "^");
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)
1150 if (seg_idx < vp->n_segs)
1151 y = vp->segs[seg_idx].points[0].y;
1152 /* else we're done */
1156 asi = active_segs[0];
1157 if (n_ips[asi] == 1)
1158 y = vp->segs[asi].points[cursor[asi] + 1].y;
1161 for (i = 1; i < n_active_segs; i++)
1163 asi = active_segs[i];
1164 if (n_ips[asi] == 1)
1165 y_curs = vp->segs[asi].points[cursor[asi] + 1].y;
1167 y_curs = ips[asi][1].y;
1171 if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
1172 y = vp->segs[seg_idx].points[0].y;
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++)
1181 asi = active_segs[i];
1182 if (n_ips[asi] == 1)
1183 p_curs = vp->segs[asi].points[cursor[asi] + 1];
1185 p_curs = ips[asi][1];
1188 svp_add_point (new_vp, n_points_max,
1189 p_curs, seg_map, active_segs, n_active_segs, i);
1192 for (j = 0; j < n_ips[asi]; j++)
1193 ips[asi][j] = ips[asi][j + 1];
1195 if (n_ips[asi] == 0)
1197 ips[asi][0] = p_curs;
1202 if (first_share < 0 || p_curs.x != share_x)
1204 /* this is where crossings are detected, and if
1205 found, the active segments switched around. */
1207 fix_crossing (first_share, i,
1208 active_segs, n_active_segs,
1209 cursor, ips, n_ips, n_ips_max, vp, seg_map,
1211 &n_segs_max, &n_points_max);
1217 if (cursor[asi] < vp->segs[asi].n_points - 1)
1221 intersect_neighbors (i, active_segs,
1222 n_ips, n_ips_max, ips,
1225 if (i + 1 < n_active_segs)
1226 intersect_neighbors (i + 1, active_segs,
1227 n_ips, n_ips_max, ips,
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,
1238 &n_segs_max, &n_points_max);
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,
1248 &n_segs_max, &n_points_max);
1255 /* not necessary to sort, new segments only get added at y, which
1256 increases monotonically */
1258 qsort (&new_vp->segs, new_vp->n_segs, sizeof (svp_seg), svp_seg_compare);
1261 for (k = 0; k < new_vp->n_segs - 1; k++)
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);
1277 art_free (n_points_max);
1279 art_free (n_ips_max);
1283 art_free (active_segs);
1290 /* Rewind a svp satisfying the nocross invariant.
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.
1297 Status of this routine:
1299 Basic correctness: Was ok in gfonted. However, this code does not
1300 yet compute bboxes for the resulting svp segs.
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.
1308 Precision: good, except that no attempt is made to remove "hair".
1309 Doing random perturbation just makes matters worse.
1313 * art_svp_rewind_uncrossed: Rewind an svp satisfying the nocross invariant.
1314 * @vp: The original svp.
1315 * @rule: The winding rule.
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.
1321 * With @rule == ART_WIND_RULE_NONZERO, the resulting svp has a
1322 * winding number of 1 where @vp has a nonzero winding number.
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.
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.
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.
1337 * This routine needs to be redone from scratch with numerical robustness
1338 * in mind. I'm working on it.
1340 * Return value: The new svp.
1343 art_svp_rewind_uncrossed (ArtSVP *vp, ArtWindRule rule)
1365 new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) +
1366 (n_segs_max - 1) * sizeof(ArtSVPSeg));
1369 if (vp->n_segs == 0)
1372 winding = art_new (int, vp->n_segs);
1374 active_segs = art_new (int, vp->n_segs);
1375 cursor = art_new (int, vp->n_segs);
1379 y = vp->segs[0].points[0].y;
1380 while (seg_idx < vp->n_segs || n_active_segs > 0)
1383 printf ("y = %g\n", y);
1385 /* delete segments ending at y from active list */
1386 for (i = 0; i < n_active_segs; i++)
1388 asi = active_segs[i];
1389 if (vp->segs[asi].n_points - 1 == cursor[asi] &&
1390 vp->segs[asi].points[cursor[asi]].y == y)
1393 printf ("deleting %d\n", asi);
1396 for (j = i; j < n_active_segs; j++)
1397 active_segs[j] = active_segs[j + 1];
1402 /* insert new segments into the active list */
1403 while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y)
1406 printf ("inserting %d\n", seg_idx);
1408 cursor[seg_idx] = 0;
1409 for (i = 0; i < n_active_segs; i++)
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)
1419 /* Determine winding number for this segment */
1422 else if (vp->segs[active_segs[i - 1]].dir)
1423 left_wind = winding[active_segs[i - 1]];
1425 left_wind = winding[active_segs[i - 1]] - 1;
1427 if (vp->segs[seg_idx].dir)
1428 wind = left_wind + 1;
1432 winding[seg_idx] = wind;
1436 case ART_WIND_RULE_NONZERO:
1437 keep = (wind == 1 || wind == 0);
1438 invert = (wind == 0);
1440 case ART_WIND_RULE_INTERSECT:
1444 case ART_WIND_RULE_ODDEVEN:
1446 invert = !(wind & 1);
1448 case ART_WIND_RULE_POSITIVE:
1460 ArtPoint *points, *new_points;
1465 printf ("keeping segment %d\n", seg_idx);
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,
1474 n_points, new_dir, new_points,
1475 &vp->segs[seg_idx].bbox);
1479 for (j = i; j < n_active_segs; j++)
1481 tmp2 = active_segs[j];
1482 active_segs[j] = tmp1;
1485 active_segs[n_active_segs] = tmp1;
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++)
1495 asi = active_segs[i];
1496 printf ("%d:%d (%g, %g) - (%g, %g) %s %d\n", 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" : "^",
1507 /* advance y to the next event */
1508 if (n_active_segs == 0)
1510 if (seg_idx < vp->n_segs)
1511 y = vp->segs[seg_idx].points[0].y;
1512 /* else we're done */
1516 asi = active_segs[0];
1517 y = vp->segs[asi].points[cursor[asi] + 1].y;
1518 for (i = 1; i < n_active_segs; i++)
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;
1524 if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y)
1525 y = vp->segs[seg_idx].points[0].y;
1528 /* advance cursors to reach new y */
1529 for (i = 0; i < n_active_segs; i++)
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)
1541 art_free (active_segs);