X-Git-Url: http://git.asbjorn.biz/?a=blobdiff_plain;f=lib%2Fart%2Fart_svp_wind.c;fp=lib%2Fart%2Fart_svp_wind.c;h=a12b1c7763a770c5c94edf9babe3ac6463eedaff;hb=418ef05b5d0b74e6f19200ffa40bcc3afccb5029;hp=0000000000000000000000000000000000000000;hpb=c1f4a906117e98e2f4ba59abc1ce153e8129bc19;p=swftools.git diff --git a/lib/art/art_svp_wind.c b/lib/art/art_svp_wind.c new file mode 100644 index 0000000..a12b1c7 --- /dev/null +++ b/lib/art/art_svp_wind.c @@ -0,0 +1,1545 @@ +/* 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; +}