/* Libart_LGPL - library of basic graphic primitives * Copyright (C) 1998-2000 Raph Levien * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Library General Public License for more details. * * You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ /* Primitive intersection and winding number operations on sorted vector paths. These routines are internal to libart, used to construct operations like intersection, union, and difference. */ #include "config.h" #include "art_svp_wind.h" #include /* for printf of debugging info */ #include /* for memcpy */ #include #include "art_misc.h" #include "art_rect.h" #include "art_svp.h" #define noVERBOSE #define PT_EQ(p1,p2) ((p1).x == (p2).x && (p1).y == (p2).y) #define PT_CLOSE(p1,p2) (fabs ((p1).x - (p2).x) < 1e-6 && fabs ((p1).y - (p2).y) < 1e-6) /* return nonzero and set *p to the intersection point if the lines z0-z1 and z2-z3 intersect each other. */ static int intersect_lines (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3, ArtPoint *p) { double a01, b01, c01; double a23, b23, c23; double d0, d1, d2, d3; double det; /* if the vectors share an endpoint, they don't intersect */ if (PT_EQ (z0, z2) || PT_EQ (z0, z3) || PT_EQ (z1, z2) || PT_EQ (z1, z3)) return 0; #if 0 if (PT_CLOSE (z0, z2) || PT_CLOSE (z0, z3) || PT_CLOSE (z1, z2) || PT_CLOSE (z1, z3)) return 0; #endif /* find line equations ax + by + c = 0 */ a01 = z0.y - z1.y; b01 = z1.x - z0.x; c01 = -(z0.x * a01 + z0.y * b01); /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y) = (z1.x * z0.y - z1.y * z0.x) */ d2 = a01 * z2.x + b01 * z2.y + c01; d3 = a01 * z3.x + b01 * z3.y + c01; if ((d2 > 0) == (d3 > 0)) return 0; a23 = z2.y - z3.y; b23 = z3.x - z2.x; c23 = -(z2.x * a23 + z2.y * b23); d0 = a23 * z0.x + b23 * z0.y + c23; d1 = a23 * z1.x + b23 * z1.y + c23; if ((d0 > 0) == (d1 > 0)) return 0; /* now we definitely know that the lines intersect */ /* solve the two linear equations ax + by + c = 0 */ det = 1.0 / (a01 * b23 - a23 * b01); p->x = det * (c23 * b01 - c01 * b23); p->y = det * (c01 * a23 - c23 * a01); return 1; } #define EPSILON 1e-6 static double trap_epsilon (double v) { const double epsilon = EPSILON; if (v < epsilon && v > -epsilon) return 0; else return v; } /* Determine the order of line segments z0-z1 and z2-z3. Return +1 if z2-z3 lies entirely to the right of z0-z1, -1 if entirely to the left, or 0 if overlap. The case analysis in this function is quite ugly. The fact that it's almost 200 lines long is ridiculous. Ok, so here's the plan to cut it down: First, do a bounding line comparison on the x coordinates. This is pretty much the common case, and should go quickly. It also takes care of the case where both lines are horizontal. Then, do d0 and d1 computation, but only if a23 is nonzero. Finally, do d2 and d3 computation, but only if a01 is nonzero. Fall through to returning 0 (this will happen when both lines are horizontal and they overlap). */ static int x_order (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3) { double a01, b01, c01; double a23, b23, c23; double d0, d1, d2, d3; if (z0.y == z1.y) { if (z2.y == z3.y) { double x01min, x01max; double x23min, x23max; if (z0.x > z1.x) { x01min = z1.x; x01max = z0.x; } else { x01min = z0.x; x01max = z1.x; } if (z2.x > z3.x) { x23min = z3.x; x23max = z2.x; } else { x23min = z2.x; x23max = z3.x; } if (x23min >= x01max) return 1; else if (x01min >= x23max) return -1; else return 0; } else { /* z0-z1 is horizontal, z2-z3 isn't */ a23 = z2.y - z3.y; b23 = z3.x - z2.x; c23 = -(z2.x * a23 + z2.y * b23); if (z3.y < z2.y) { a23 = -a23; b23 = -b23; c23 = -c23; } d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23); d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23); if (d0 > 0) { if (d1 >= 0) return 1; else return 0; } else if (d0 == 0) { if (d1 > 0) return 1; else if (d1 < 0) return -1; else printf ("case 1 degenerate\n"); return 0; } else /* d0 < 0 */ { if (d1 <= 0) return -1; else return 0; } } } else if (z2.y == z3.y) { /* z2-z3 is horizontal, z0-z1 isn't */ a01 = z0.y - z1.y; b01 = z1.x - z0.x; c01 = -(z0.x * a01 + z0.y * b01); /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y) = (z1.x * z0.y - z1.y * z0.x) */ if (z1.y < z0.y) { a01 = -a01; b01 = -b01; c01 = -c01; } d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01); d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01); if (d2 > 0) { if (d3 >= 0) return -1; else return 0; } else if (d2 == 0) { if (d3 > 0) return -1; else if (d3 < 0) return 1; else printf ("case 2 degenerate\n"); return 0; } else /* d2 < 0 */ { if (d3 <= 0) return 1; else return 0; } } /* find line equations ax + by + c = 0 */ a01 = z0.y - z1.y; b01 = z1.x - z0.x; c01 = -(z0.x * a01 + z0.y * b01); /* = -((z0.y - z1.y) * z0.x + (z1.x - z0.x) * z0.y) = -(z1.x * z0.y - z1.y * z0.x) */ if (a01 > 0) { a01 = -a01; b01 = -b01; c01 = -c01; } /* so now, (a01, b01) points to the left, thus a01 * x + b01 * y + c01 is negative if the point lies to the right of the line */ d2 = trap_epsilon (a01 * z2.x + b01 * z2.y + c01); d3 = trap_epsilon (a01 * z3.x + b01 * z3.y + c01); if (d2 > 0) { if (d3 >= 0) return -1; } else if (d2 == 0) { if (d3 > 0) return -1; else if (d3 < 0) return 1; else fprintf (stderr, "colinear!\n"); } else /* d2 < 0 */ { if (d3 <= 0) return 1; } a23 = z2.y - z3.y; b23 = z3.x - z2.x; c23 = -(z2.x * a23 + z2.y * b23); if (a23 > 0) { a23 = -a23; b23 = -b23; c23 = -c23; } d0 = trap_epsilon (a23 * z0.x + b23 * z0.y + c23); d1 = trap_epsilon (a23 * z1.x + b23 * z1.y + c23); if (d0 > 0) { if (d1 >= 0) return 1; } else if (d0 == 0) { if (d1 > 0) return 1; else if (d1 < 0) return -1; else fprintf (stderr, "colinear!\n"); } else /* d0 < 0 */ { if (d1 <= 0) return -1; } return 0; } /* similar to x_order, but to determine whether point z0 + epsilon lies to the left of the line z2-z3 or to the right */ static int x_order_2 (ArtPoint z0, ArtPoint z1, ArtPoint z2, ArtPoint z3) { double a23, b23, c23; double d0, d1; a23 = z2.y - z3.y; b23 = z3.x - z2.x; c23 = -(z2.x * a23 + z2.y * b23); if (a23 > 0) { a23 = -a23; b23 = -b23; c23 = -c23; } d0 = a23 * z0.x + b23 * z0.y + c23; if (d0 > EPSILON) return -1; else if (d0 < -EPSILON) return 1; d1 = a23 * z1.x + b23 * z1.y + c23; if (d1 > EPSILON) return -1; else if (d1 < -EPSILON) return 1; if (z0.x == z1.x && z1.x == z2.x && z2.x == z3.x) { art_dprint ("x_order_2: colinear and horizontally aligned!\n"); return 0; } if (z0.x <= z2.x && z1.x <= z2.x && z0.x <= z3.x && z1.x <= z3.x) return -1; if (z0.x >= z2.x && z1.x >= z2.x && z0.x >= z3.x && z1.x >= z3.x) return 1; fprintf (stderr, "x_order_2: colinear!\n"); return 0; } #ifdef DEAD_CODE /* Traverse the vector path, keeping it in x-sorted order. This routine doesn't actually do anything - it's just here for explanatory purposes. */ void traverse (ArtSVP *vp) { int *active_segs; int n_active_segs; int *cursor; int seg_idx; double y; int tmp1, tmp2; int asi; int i, j; active_segs = art_new (int, vp->n_segs); cursor = art_new (int, vp->n_segs); n_active_segs = 0; seg_idx = 0; y = vp->segs[0].points[0].y; while (seg_idx < vp->n_segs || n_active_segs > 0) { printf ("y = %g\n", y); /* delete segments ending at y from active list */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; if (vp->segs[asi].n_points - 1 == cursor[asi] && vp->segs[asi].points[cursor[asi]].y == y) { printf ("deleting %d\n", asi); n_active_segs--; for (j = i; j < n_active_segs; j++) active_segs[j] = active_segs[j + 1]; i--; } } /* insert new segments into the active list */ while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y) { cursor[seg_idx] = 0; printf ("inserting %d\n", seg_idx); for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; if (x_order (vp->segs[asi].points[cursor[asi]], vp->segs[asi].points[cursor[asi] + 1], vp->segs[seg_idx].points[0], vp->segs[seg_idx].points[1]) == -1) break; } tmp1 = seg_idx; for (j = i; j < n_active_segs; j++) { tmp2 = active_segs[j]; active_segs[j] = tmp1; tmp1 = tmp2; } active_segs[n_active_segs] = tmp1; n_active_segs++; seg_idx++; } /* all active segs cross the y scanline (considering segs to be closed on top and open on bottom) */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; printf ("%d (%g, %g) - (%g, %g) %s\n", asi, vp->segs[asi].points[cursor[asi]].x, vp->segs[asi].points[cursor[asi]].y, vp->segs[asi].points[cursor[asi] + 1].x, vp->segs[asi].points[cursor[asi] + 1].y, vp->segs[asi].dir ? "v" : "^"); } /* advance y to the next event */ if (n_active_segs == 0) { if (seg_idx < vp->n_segs) y = vp->segs[seg_idx].points[0].y; /* else we're done */ } else { asi = active_segs[0]; y = vp->segs[asi].points[cursor[asi] + 1].y; for (i = 1; i < n_active_segs; i++) { asi = active_segs[i]; if (y > vp->segs[asi].points[cursor[asi] + 1].y) y = vp->segs[asi].points[cursor[asi] + 1].y; } if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y) y = vp->segs[seg_idx].points[0].y; } /* advance cursors to reach new y */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; while (cursor[asi] < vp->segs[asi].n_points - 1 && y >= vp->segs[asi].points[cursor[asi] + 1].y) cursor[asi]++; } printf ("\n"); } art_free (cursor); art_free (active_segs); } #endif /* I believe that the loop will always break with i=1. I think I'll want to change this from a simple sorted list to a modified stack. ips[*][0] will get its own data structure, and ips[*] will in general only be allocated if there is an intersection. Finally, the segment can be traced through the initial point (formerly ips[*][0]), backwards through the stack, and finally to cursor + 1. This change should cut down on allocation bandwidth, and also eliminate the iteration through n_ipl below. */ static void insert_ip (int seg_i, int *n_ips, int *n_ips_max, ArtPoint **ips, ArtPoint ip) { int i; ArtPoint tmp1, tmp2; int n_ipl; ArtPoint *ipl; n_ipl = n_ips[seg_i]++; if (n_ipl == n_ips_max[seg_i]) art_expand (ips[seg_i], ArtPoint, n_ips_max[seg_i]); ipl = ips[seg_i]; for (i = 1; i < n_ipl; i++) if (ipl[i].y > ip.y) break; tmp1 = ip; for (; i <= n_ipl; i++) { tmp2 = ipl[i]; ipl[i] = tmp1; tmp1 = tmp2; } } /* test active segment (i - 1) against i for intersection, if so, add intersection point to both ips lists. */ static void intersect_neighbors (int i, int *active_segs, int *n_ips, int *n_ips_max, ArtPoint **ips, int *cursor, ArtSVP *vp) { ArtPoint z0, z1, z2, z3; int asi01, asi23; ArtPoint ip; asi01 = active_segs[i - 1]; z0 = ips[asi01][0]; if (n_ips[asi01] == 1) z1 = vp->segs[asi01].points[cursor[asi01] + 1]; else z1 = ips[asi01][1]; asi23 = active_segs[i]; z2 = ips[asi23][0]; if (n_ips[asi23] == 1) z3 = vp->segs[asi23].points[cursor[asi23] + 1]; else z3 = ips[asi23][1]; if (intersect_lines (z0, z1, z2, z3, &ip)) { #ifdef VERBOSE printf ("new intersection point: (%g, %g)\n", ip.x, ip.y); #endif insert_ip (asi01, n_ips, n_ips_max, ips, ip); insert_ip (asi23, n_ips, n_ips_max, ips, ip); } } /* Add a new point to a segment in the svp. Here, we also check to make sure that the segments satisfy nocross. However, this is only valuable for debugging, and could possibly be removed. */ static void svp_add_point (ArtSVP *svp, int *n_points_max, ArtPoint p, int *seg_map, int *active_segs, int n_active_segs, int i) { int asi, asi_left, asi_right; int n_points, n_points_left, n_points_right; ArtSVPSeg *seg; asi = seg_map[active_segs[i]]; seg = &svp->segs[asi]; n_points = seg->n_points; /* find out whether neighboring segments share a point */ if (i > 0) { asi_left = seg_map[active_segs[i - 1]]; n_points_left = svp->segs[asi_left].n_points; if (n_points_left > 1 && PT_EQ (svp->segs[asi_left].points[n_points_left - 2], svp->segs[asi].points[n_points - 1])) { /* ok, new vector shares a top point with segment to the left - now, check that it satisfies ordering invariant */ if (x_order (svp->segs[asi_left].points[n_points_left - 2], svp->segs[asi_left].points[n_points_left - 1], svp->segs[asi].points[n_points - 1], p) < 1) { #ifdef VERBOSE printf ("svp_add_point: cross on left!\n"); #endif } } } if (i + 1 < n_active_segs) { asi_right = seg_map[active_segs[i + 1]]; n_points_right = svp->segs[asi_right].n_points; if (n_points_right > 1 && PT_EQ (svp->segs[asi_right].points[n_points_right - 2], svp->segs[asi].points[n_points - 1])) { /* ok, new vector shares a top point with segment to the right - now, check that it satisfies ordering invariant */ if (x_order (svp->segs[asi_right].points[n_points_right - 2], svp->segs[asi_right].points[n_points_right - 1], svp->segs[asi].points[n_points - 1], p) > -1) { #ifdef VERBOSE printf ("svp_add_point: cross on right!\n"); #endif } } } if (n_points_max[asi] == n_points) art_expand (seg->points, ArtPoint, n_points_max[asi]); seg->points[n_points] = p; if (p.x < seg->bbox.x0) seg->bbox.x0 = p.x; else if (p.x > seg->bbox.x1) seg->bbox.x1 = p.x; seg->bbox.y1 = p.y; seg->n_points++; } #if 0 /* find where the segment (currently at i) is supposed to go, and return the target index - if equal to i, then there is no crossing problem. "Where it is supposed to go" is defined as following: Delete element i, re-insert at position target (bumping everything target and greater to the right). */ static int find_crossing (int i, int *active_segs, int n_active_segs, int *cursor, ArtPoint **ips, int *n_ips, ArtSVP *vp) { int asi, asi_left, asi_right; ArtPoint p0, p1; ArtPoint p0l, p1l; ArtPoint p0r, p1r; int target; asi = active_segs[i]; p0 = ips[asi][0]; if (n_ips[asi] == 1) p1 = vp->segs[asi].points[cursor[asi] + 1]; else p1 = ips[asi][1]; for (target = i; target > 0; target--) { asi_left = active_segs[target - 1]; p0l = ips[asi_left][0]; if (n_ips[asi_left] == 1) p1l = vp->segs[asi_left].points[cursor[asi_left] + 1]; else p1l = ips[asi_left][1]; if (!PT_EQ (p0, p0l)) break; #ifdef VERBOSE printf ("point matches on left (%g, %g) - (%g, %g) x (%g, %g) - (%g, %g)!\n", p0l.x, p0l.y, p1l.x, p1l.y, p0.x, p0.y, p1.x, p1.y); #endif if (x_order (p0l, p1l, p0, p1) == 1) break; #ifdef VERBOSE printf ("scanning to the left (i=%d, target=%d)\n", i, target); #endif } if (target < i) return target; for (; target < n_active_segs - 1; target++) { asi_right = active_segs[target + 1]; p0r = ips[asi_right][0]; if (n_ips[asi_right] == 1) p1r = vp->segs[asi_right].points[cursor[asi_right] + 1]; else p1r = ips[asi_right][1]; if (!PT_EQ (p0, p0r)) break; #ifdef VERBOSE printf ("point matches on left (%g, %g) - (%g, %g) x (%g, %g) - (%g, %g)!\n", p0.x, p0.y, p1.x, p1.y, p0r.x, p0r.y, p1r.x, p1r.y); #endif if (x_order (p0r, p1r, p0, p1) == 1) break; #ifdef VERBOSE printf ("scanning to the right (i=%d, target=%d)\n", i, target); #endif } return target; } #endif /* This routine handles the case where the segment changes its position in the active segment list. Generally, this will happen when the segment (defined by i and cursor) shares a top point with a neighbor, but breaks the ordering invariant. Essentially, this routine sorts the lines [start..end), all of which share a top point. This is implemented as your basic insertion sort. This routine takes care of intersecting the appropriate neighbors, as well. A first argument of -1 immediately returns, which helps reduce special casing in the main unwind routine. */ static void fix_crossing (int start, int end, int *active_segs, int n_active_segs, int *cursor, ArtPoint **ips, int *n_ips, int *n_ips_max, ArtSVP *vp, int *seg_map, ArtSVP **p_new_vp, int *pn_segs_max, int **pn_points_max) { int i, j; int target; int asi, asj; ArtPoint p0i, p1i; ArtPoint p0j, p1j; int swap = 0; #ifdef VERBOSE int k; #endif ArtPoint *pts; #ifdef VERBOSE printf ("fix_crossing: [%d..%d)", start, end); for (k = 0; k < n_active_segs; k++) printf (" %d", active_segs[k]); printf ("\n"); #endif if (start == -1) return; for (i = start + 1; i < end; i++) { asi = active_segs[i]; if (cursor[asi] < vp->segs[asi].n_points - 1) { p0i = ips[asi][0]; if (n_ips[asi] == 1) p1i = vp->segs[asi].points[cursor[asi] + 1]; else p1i = ips[asi][1]; for (j = i - 1; j >= start; j--) { asj = active_segs[j]; if (cursor[asj] < vp->segs[asj].n_points - 1) { p0j = ips[asj][0]; if (n_ips[asj] == 1) p1j = vp->segs[asj].points[cursor[asj] + 1]; else p1j = ips[asj][1]; /* we _hope_ p0i = p0j */ if (x_order_2 (p0j, p1j, p0i, p1i) == -1) break; } } target = j + 1; /* target is where active_seg[i] _should_ be in active_segs */ if (target != i) { swap = 1; #ifdef VERBOSE printf ("fix_crossing: at %i should be %i\n", i, target); #endif /* let's close off all relevant segments */ for (j = i; j >= target; j--) { asi = active_segs[j]; /* First conjunct: this isn't the last point in the original segment. Second conjunct: this isn't the first point in the new segment (i.e. already broken). */ if (cursor[asi] < vp->segs[asi].n_points - 1 && (*p_new_vp)->segs[seg_map[asi]].n_points != 1) { int seg_num; /* so break here */ #ifdef VERBOSE printf ("closing off %d\n", j); #endif pts = art_new (ArtPoint, 16); pts[0] = ips[asi][0]; seg_num = art_svp_add_segment (p_new_vp, pn_segs_max, pn_points_max, 1, vp->segs[asi].dir, pts, NULL); (*pn_points_max)[seg_num] = 16; seg_map[asi] = seg_num; } } /* now fix the ordering in active_segs */ asi = active_segs[i]; for (j = i; j > target; j--) active_segs[j] = active_segs[j - 1]; active_segs[j] = asi; } } } if (swap && start > 0) { int as_start; as_start = active_segs[start]; if (cursor[as_start] < vp->segs[as_start].n_points) { #ifdef VERBOSE printf ("checking intersection of %d, %d\n", start - 1, start); #endif intersect_neighbors (start, active_segs, n_ips, n_ips_max, ips, cursor, vp); } } if (swap && end < n_active_segs) { int as_end; as_end = active_segs[end - 1]; if (cursor[as_end] < vp->segs[as_end].n_points) { #ifdef VERBOSE printf ("checking intersection of %d, %d\n", end - 1, end); #endif intersect_neighbors (end, active_segs, n_ips, n_ips_max, ips, cursor, vp); } } if (swap) { #ifdef VERBOSE printf ("fix_crossing return: [%d..%d)", start, end); for (k = 0; k < n_active_segs; k++) printf (" %d", active_segs[k]); printf ("\n"); #endif } } /* Return a new sorted vector that covers the same area as the argument, but which satisfies the nocross invariant. Basically, this routine works by finding the intersection points, and cutting the segments at those points. Status of this routine: Basic correctness: Seems ok. Numerical stability: known problems in the case of points falling on lines, and colinear lines. For actual use, randomly perturbing the vertices is currently recommended. Speed: pretty good, although a more efficient priority queue, as well as bbox culling of potential intersections, are two optimizations that could help. Precision: pretty good, although the numerical stability problems make this routine unsuitable for precise calculations of differences. */ /* Here is a more detailed description of the algorithm. It follows roughly the structure of traverse (above), but is obviously quite a bit more complex. Here are a few important data structures: A new sorted vector path (new_svp). For each (active) segment in the original, a list of intersection points. Of course, the original being traversed. The following invariants hold (in addition to the invariants of the traverse procedure). The new sorted vector path lies entirely above the y scan line. The new sorted vector path keeps the nocross invariant. For each active segment, the y scan line crosses the line from the first to the second of the intersection points (where the second point is cursor + 1 if there is only one intersection point). The list of intersection points + the (cursor + 1) point is kept in nondecreasing y order. Of the active segments, none of the lines from first to second intersection point cross the 1st ip..2nd ip line of the left or right neighbor. (However, such a line may cross further intersection points of the neighbors, or segments past the immediate neighbors). Of the active segments, all lines from 1st ip..2nd ip are in strictly increasing x_order (this is very similar to the invariant of the traverse procedure, but is explicitly stated here in terms of ips). (this basically says that nocross holds on the active segments) The combination of the new sorted vector path, the path through all the intersection points to cursor + 1, and [cursor + 1, n_points) covers the same area as the argument. Another important data structure is mapping from original segment number to new segment number. The algorithm is perhaps best understood as advancing the cursors while maintaining these invariants. Here's roughly how it's done. When deleting segments from the active list, those segments are added to the new sorted vector path. In addition, the neighbors may intersect each other, so they are intersection tested (see below). When inserting new segments, they are intersection tested against their neighbors. The top point of the segment becomes the first intersection point. Advancing the cursor is just a bit different from the traverse routine, as the cursor may advance through the intersection points as well. Only when there is a single intersection point in the list does the cursor advance in the original segment. In either case, the new vector is intersection tested against both neighbors. It also causes the vector over which the cursor is advancing to be added to the new svp. Two steps need further clarification: Intersection testing: the 1st ip..2nd ip lines of the neighbors are tested to see if they cross (using intersect_lines). If so, then the intersection point is added to the ip list of both segments, maintaining the invariant that the list of intersection points is nondecreasing in y). Adding vector to new svp: if the new vector shares a top x coordinate with another vector, then it is checked to see whether it is in order. If not, then both segments are "broken," and then restarted. Note: in the case when both segments are in the same order, they may simply be swapped without breaking. For the time being, I'm going to put some of these operations into subroutines. If it turns out to be a performance problem, I could try to reorganize the traverse procedure so that each is only called once, and inline them. But if it's not a performance problem, I'll just keep it this way, because it will probably help to make the code clearer, and I believe this code could use all the clarity it can get. */ /** * art_svp_uncross: Resolve self-intersections of an svp. * @vp: The original svp. * * Finds all the intersections within @vp, and constructs a new svp * with new points added at these intersections. * * This routine needs to be redone from scratch with numerical robustness * in mind. I'm working on it. * * Return value: The new svp. **/ ArtSVP * art_svp_uncross (ArtSVP *vp) { int *active_segs; int n_active_segs; int *cursor; int seg_idx; double y; int tmp1, tmp2; int asi; int i, j; /* new data structures */ /* intersection points; invariant: *ips[i] is only allocated if i is active */ int *n_ips, *n_ips_max; ArtPoint **ips; /* new sorted vector path */ int n_segs_max, seg_num; ArtSVP *new_vp; int *n_points_max; /* mapping from argument to new segment numbers - again, only valid if active */ int *seg_map; double y_curs; ArtPoint p_curs; int first_share; double share_x; ArtPoint *pts; n_segs_max = 16; new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) + (n_segs_max - 1) * sizeof(ArtSVPSeg)); new_vp->n_segs = 0; if (vp->n_segs == 0) return new_vp; active_segs = art_new (int, vp->n_segs); cursor = art_new (int, vp->n_segs); seg_map = art_new (int, vp->n_segs); n_ips = art_new (int, vp->n_segs); n_ips_max = art_new (int, vp->n_segs); ips = art_new (ArtPoint *, vp->n_segs); n_points_max = art_new (int, n_segs_max); n_active_segs = 0; seg_idx = 0; y = vp->segs[0].points[0].y; while (seg_idx < vp->n_segs || n_active_segs > 0) { #ifdef VERBOSE printf ("y = %g\n", y); #endif /* maybe move deletions to end of loop (to avoid so much special casing on the end of a segment)? */ /* delete segments ending at y from active list */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; if (vp->segs[asi].n_points - 1 == cursor[asi] && vp->segs[asi].points[cursor[asi]].y == y) { do { #ifdef VERBOSE printf ("deleting %d\n", asi); #endif art_free (ips[asi]); n_active_segs--; for (j = i; j < n_active_segs; j++) active_segs[j] = active_segs[j + 1]; if (i < n_active_segs) asi = active_segs[i]; else break; } while (vp->segs[asi].n_points - 1 == cursor[asi] && vp->segs[asi].points[cursor[asi]].y == y); /* test intersection of neighbors */ if (i > 0 && i < n_active_segs) intersect_neighbors (i, active_segs, n_ips, n_ips_max, ips, cursor, vp); i--; } } /* insert new segments into the active list */ while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y) { #ifdef VERBOSE printf ("inserting %d\n", seg_idx); #endif cursor[seg_idx] = 0; for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; if (x_order_2 (vp->segs[seg_idx].points[0], vp->segs[seg_idx].points[1], vp->segs[asi].points[cursor[asi]], vp->segs[asi].points[cursor[asi] + 1]) == -1) break; } /* Create and initialize the intersection points data structure */ n_ips[seg_idx] = 1; n_ips_max[seg_idx] = 2; ips[seg_idx] = art_new (ArtPoint, n_ips_max[seg_idx]); ips[seg_idx][0] = vp->segs[seg_idx].points[0]; /* Start a new segment in the new vector path */ pts = art_new (ArtPoint, 16); pts[0] = vp->segs[seg_idx].points[0]; seg_num = art_svp_add_segment (&new_vp, &n_segs_max, &n_points_max, 1, vp->segs[seg_idx].dir, pts, NULL); n_points_max[seg_num] = 16; seg_map[seg_idx] = seg_num; tmp1 = seg_idx; for (j = i; j < n_active_segs; j++) { tmp2 = active_segs[j]; active_segs[j] = tmp1; tmp1 = tmp2; } active_segs[n_active_segs] = tmp1; n_active_segs++; if (i > 0) intersect_neighbors (i, active_segs, n_ips, n_ips_max, ips, cursor, vp); if (i + 1 < n_active_segs) intersect_neighbors (i + 1, active_segs, n_ips, n_ips_max, ips, cursor, vp); seg_idx++; } /* all active segs cross the y scanline (considering segs to be closed on top and open on bottom) */ #ifdef VERBOSE for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; printf ("%d ", asi); for (j = 0; j < n_ips[asi]; j++) printf ("(%g, %g) - ", ips[asi][j].x, ips[asi][j].y); printf ("(%g, %g) %s\n", vp->segs[asi].points[cursor[asi] + 1].x, vp->segs[asi].points[cursor[asi] + 1].y, vp->segs[asi].dir ? "v" : "^"); } #endif /* advance y to the next event Note: this is quadratic. We'd probably get decent constant factor speed improvement by caching the y_curs values. */ if (n_active_segs == 0) { if (seg_idx < vp->n_segs) y = vp->segs[seg_idx].points[0].y; /* else we're done */ } else { asi = active_segs[0]; if (n_ips[asi] == 1) y = vp->segs[asi].points[cursor[asi] + 1].y; else y = ips[asi][1].y; for (i = 1; i < n_active_segs; i++) { asi = active_segs[i]; if (n_ips[asi] == 1) y_curs = vp->segs[asi].points[cursor[asi] + 1].y; else y_curs = ips[asi][1].y; if (y > y_curs) y = y_curs; } if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y) y = vp->segs[seg_idx].points[0].y; } first_share = -1; share_x = 0; /* to avoid gcc warning, although share_x is never used when first_share is -1 */ /* advance cursors to reach new y */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; if (n_ips[asi] == 1) p_curs = vp->segs[asi].points[cursor[asi] + 1]; else p_curs = ips[asi][1]; if (p_curs.y == y) { svp_add_point (new_vp, n_points_max, p_curs, seg_map, active_segs, n_active_segs, i); n_ips[asi]--; for (j = 0; j < n_ips[asi]; j++) ips[asi][j] = ips[asi][j + 1]; if (n_ips[asi] == 0) { ips[asi][0] = p_curs; n_ips[asi] = 1; cursor[asi]++; } if (first_share < 0 || p_curs.x != share_x) { /* this is where crossings are detected, and if found, the active segments switched around. */ fix_crossing (first_share, i, active_segs, n_active_segs, cursor, ips, n_ips, n_ips_max, vp, seg_map, &new_vp, &n_segs_max, &n_points_max); first_share = i; share_x = p_curs.x; } if (cursor[asi] < vp->segs[asi].n_points - 1) { if (i > 0) intersect_neighbors (i, active_segs, n_ips, n_ips_max, ips, cursor, vp); if (i + 1 < n_active_segs) intersect_neighbors (i + 1, active_segs, n_ips, n_ips_max, ips, cursor, vp); } } else { /* not on a cursor point */ fix_crossing (first_share, i, active_segs, n_active_segs, cursor, ips, n_ips, n_ips_max, vp, seg_map, &new_vp, &n_segs_max, &n_points_max); first_share = -1; } } /* fix crossing on last shared group */ fix_crossing (first_share, i, active_segs, n_active_segs, cursor, ips, n_ips, n_ips_max, vp, seg_map, &new_vp, &n_segs_max, &n_points_max); #ifdef VERBOSE printf ("\n"); #endif } /* not necessary to sort, new segments only get added at y, which increases monotonically */ #if 0 qsort (&new_vp->segs, new_vp->n_segs, sizeof (svp_seg), svp_seg_compare); { int k; for (k = 0; k < new_vp->n_segs - 1; k++) { printf ("(%g, %g) - (%g, %g) %s (%g, %g) - (%g, %g)\n", new_vp->segs[k].points[0].x, new_vp->segs[k].points[0].y, new_vp->segs[k].points[1].x, new_vp->segs[k].points[1].y, svp_seg_compare (&new_vp->segs[k], &new_vp->segs[k + 1]) > 1 ? ">": "<", new_vp->segs[k + 1].points[0].x, new_vp->segs[k + 1].points[0].y, new_vp->segs[k + 1].points[1].x, new_vp->segs[k + 1].points[1].y); } } #endif art_free (n_points_max); art_free (seg_map); art_free (n_ips_max); art_free (n_ips); art_free (ips); art_free (cursor); art_free (active_segs); return new_vp; } #define noVERBOSE /* Rewind a svp satisfying the nocross invariant. The winding number of a segment is defined as the winding number of the points to the left while travelling in the direction of the segment. Therefore it preincrements and postdecrements as a scan line is traversed from left to right. Status of this routine: Basic correctness: Was ok in gfonted. However, this code does not yet compute bboxes for the resulting svp segs. Numerical stability: known problems in the case of horizontal segments in polygons with any complexity. For actual use, randomly perturbing the vertices is recommended. Speed: good. Precision: good, except that no attempt is made to remove "hair". Doing random perturbation just makes matters worse. */ /** * art_svp_rewind_uncrossed: Rewind an svp satisfying the nocross invariant. * @vp: The original svp. * @rule: The winding rule. * * Creates a new svp with winding number of 0 or 1 everywhere. The @rule * argument specifies a rule for how winding numbers in the original * @vp map to the winding numbers in the result. * * With @rule == ART_WIND_RULE_NONZERO, the resulting svp has a * winding number of 1 where @vp has a nonzero winding number. * * With @rule == ART_WIND_RULE_INTERSECT, the resulting svp has a * winding number of 1 where @vp has a winding number greater than * 1. It is useful for computing intersections. * * With @rule == ART_WIND_RULE_ODDEVEN, the resulting svp has a * winding number of 1 where @vp has an odd winding number. It is * useful for implementing the even-odd winding rule of the * PostScript imaging model. * * With @rule == ART_WIND_RULE_POSITIVE, the resulting svp has a * winding number of 1 where @vp has a positive winding number. It is * useful for implementing asymmetric difference. * * This routine needs to be redone from scratch with numerical robustness * in mind. I'm working on it. * * Return value: The new svp. **/ ArtSVP * art_svp_rewind_uncrossed (ArtSVP *vp, ArtWindRule rule) { int *active_segs; int n_active_segs; int *cursor; int seg_idx; double y; int tmp1, tmp2; int asi; int i, j; ArtSVP *new_vp; int n_segs_max; int *winding; int left_wind; int wind; int keep, invert; #ifdef VERBOSE print_svp (vp); #endif n_segs_max = 16; new_vp = (ArtSVP *)art_alloc (sizeof(ArtSVP) + (n_segs_max - 1) * sizeof(ArtSVPSeg)); new_vp->n_segs = 0; if (vp->n_segs == 0) return new_vp; winding = art_new (int, vp->n_segs); active_segs = art_new (int, vp->n_segs); cursor = art_new (int, vp->n_segs); n_active_segs = 0; seg_idx = 0; y = vp->segs[0].points[0].y; while (seg_idx < vp->n_segs || n_active_segs > 0) { #ifdef VERBOSE printf ("y = %g\n", y); #endif /* delete segments ending at y from active list */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; if (vp->segs[asi].n_points - 1 == cursor[asi] && vp->segs[asi].points[cursor[asi]].y == y) { #ifdef VERBOSE printf ("deleting %d\n", asi); #endif n_active_segs--; for (j = i; j < n_active_segs; j++) active_segs[j] = active_segs[j + 1]; i--; } } /* insert new segments into the active list */ while (seg_idx < vp->n_segs && y == vp->segs[seg_idx].points[0].y) { #ifdef VERBOSE printf ("inserting %d\n", seg_idx); #endif cursor[seg_idx] = 0; for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; if (x_order_2 (vp->segs[seg_idx].points[0], vp->segs[seg_idx].points[1], vp->segs[asi].points[cursor[asi]], vp->segs[asi].points[cursor[asi] + 1]) == -1) break; } /* Determine winding number for this segment */ if (i == 0) left_wind = 0; else if (vp->segs[active_segs[i - 1]].dir) left_wind = winding[active_segs[i - 1]]; else left_wind = winding[active_segs[i - 1]] - 1; if (vp->segs[seg_idx].dir) wind = left_wind + 1; else wind = left_wind; winding[seg_idx] = wind; switch (rule) { case ART_WIND_RULE_NONZERO: keep = (wind == 1 || wind == 0); invert = (wind == 0); break; case ART_WIND_RULE_INTERSECT: keep = (wind == 2); invert = 0; break; case ART_WIND_RULE_ODDEVEN: keep = 1; invert = !(wind & 1); break; case ART_WIND_RULE_POSITIVE: keep = (wind == 1); invert = 0; break; default: keep = 0; invert = 0; break; } if (keep) { ArtPoint *points, *new_points; int n_points; int new_dir; #ifdef VERBOSE printf ("keeping segment %d\n", seg_idx); #endif n_points = vp->segs[seg_idx].n_points; points = vp->segs[seg_idx].points; new_points = art_new (ArtPoint, n_points); memcpy (new_points, points, n_points * sizeof (ArtPoint)); new_dir = vp->segs[seg_idx].dir ^ invert; art_svp_add_segment (&new_vp, &n_segs_max, NULL, n_points, new_dir, new_points, &vp->segs[seg_idx].bbox); } tmp1 = seg_idx; for (j = i; j < n_active_segs; j++) { tmp2 = active_segs[j]; active_segs[j] = tmp1; tmp1 = tmp2; } active_segs[n_active_segs] = tmp1; n_active_segs++; seg_idx++; } #ifdef VERBOSE /* all active segs cross the y scanline (considering segs to be closed on top and open on bottom) */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; printf ("%d:%d (%g, %g) - (%g, %g) %s %d\n", asi, cursor[asi], vp->segs[asi].points[cursor[asi]].x, vp->segs[asi].points[cursor[asi]].y, vp->segs[asi].points[cursor[asi] + 1].x, vp->segs[asi].points[cursor[asi] + 1].y, vp->segs[asi].dir ? "v" : "^", winding[asi]); } #endif /* advance y to the next event */ if (n_active_segs == 0) { if (seg_idx < vp->n_segs) y = vp->segs[seg_idx].points[0].y; /* else we're done */ } else { asi = active_segs[0]; y = vp->segs[asi].points[cursor[asi] + 1].y; for (i = 1; i < n_active_segs; i++) { asi = active_segs[i]; if (y > vp->segs[asi].points[cursor[asi] + 1].y) y = vp->segs[asi].points[cursor[asi] + 1].y; } if (seg_idx < vp->n_segs && y > vp->segs[seg_idx].points[0].y) y = vp->segs[seg_idx].points[0].y; } /* advance cursors to reach new y */ for (i = 0; i < n_active_segs; i++) { asi = active_segs[i]; while (cursor[asi] < vp->segs[asi].n_points - 1 && y >= vp->segs[asi].points[cursor[asi] + 1].y) cursor[asi]++; } #ifdef VERBOSE printf ("\n"); #endif } art_free (cursor); art_free (active_segs); art_free (winding); return new_vp; }