/* This file is part of the GNU plotutils package. Copyright (C) 1995, 1996, 1997, 1998, 1999, 2000, 2005, Free Software Foundation, Inc. The GNU plotutils package is free software. You may redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software foundation; either version 2, or (at your option) any later version. The GNU plotutils package 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 General Public License for more details. You should have received a copy of the GNU General Public License along with the GNU plotutils package; see the file COPYING. If not, write to the Free Software Foundation, Inc., 51 Franklin St., Fifth Floor, Boston, MA 02110-1301, USA. */ /* This file contains all of libplot's low-level code for constructing and manipulating paths, both simple and compound. It includes functions that construct polygonal and Bezier approximations to a given path. They are used by Plotters whose output format does not support all of libplot's graphics primitives. */ /* E.g., _fakearc() draws polygonal approximations to circular and elliptic quarter-arcs. Each polygonal approximation will contain 2**NUM_ARC_SUBDIVISIONS line segments. Similarly, polygonal approximations to quadratic and cubic Beziers will contain no more than 2**MAX_NUM_BEZIER2_SUBDIVISIONS and 2**MAX_NUM_BEZIER3_SUBDIVISIONS line segments. However, each bisection algorithm used for drawing a Bezier normally usually its recursion based on a relative flatness criterion (see below). */ #include "sys-defines.h" #include "extern.h" #include "g_arc.h" /* for chord table */ /* Number of times a circular arc or quarter-ellipse will be recursively subdivided. Two raised to this power is the number of line segments that the polygonalization will contain. */ /* NOTE: the maximum allowed value for NUM_ARC_SUBDIVISIONS is TABULATED_ARC_SUBDIVISIONS (i.e., the size of the chordal deviation table for a quarter-circle or quarter-ellipse, defined in g_arc.h). */ #define NUM_ARC_SUBDIVISIONS 5 /* Maximum number of times quadratic and cubic Beziers will be recursively subdivided. For Beziers we use an adaptive algorithm, in which bisection stops when a specified relative flatness has been reached. But these parameters provide a hard cutoff, which overrides the relative flatness end condition. */ #define MAX_NUM_BEZIER2_SUBDIVISIONS 6 #define MAX_NUM_BEZIER3_SUBDIVISIONS 7 /* The relative flatness parameters. */ #define REL_QUAD_FLATNESS 5e-4 #define REL_CUBIC_FLATNESS 5e-4 #define DATAPOINTS_BUFSIZ PL_MAX_UNFILLED_PATH_LENGTH #define DIST(p0,p1) (sqrt( ((p0).x - (p1).x)*((p0).x - (p1).x) \ + ((p0).y - (p1).y)*((p0).y - (p1).y))) #define MIDWAY(x0, x1) (0.5 * ((x0) + (x1))) /* forward references */ static void _prepare_chord_table (double sagitta, double custom_chord_table[TABULATED_ARC_SUBDIVISIONS]); static void _fakearc (plPath *path, plPoint p0, plPoint p1, int arc_type, const double *custom_chord_table, const double m[4]); /* ctor for plPath class; constructs an empty plPath, with type set to PATH_SEGMENT_LIST (default type) */ plPath * _new_plPath (void) { plPath *path; path = (plPath *)_pl_xmalloc (sizeof (plPath)); path->type = PATH_SEGMENT_LIST; path->segments = (plPathSegment *)NULL; path->segments_len = 0; /* number of slots allocated */ path->num_segments = 0; /* number of slots occupied */ path->primitive = false; path->llx = DBL_MAX; path->lly = DBL_MAX; path->urx = -(DBL_MAX); path->ury = -(DBL_MAX); return path; } /* dtor for plPath class */ void _delete_plPath (plPath *path) { if (path == (plPath *)NULL) return; if (path->type == PATH_SEGMENT_LIST && path->segments_len > 0) /* number of slots allocated */ free (path->segments); free (path); } /* reset function for plPath class */ void _reset_plPath (plPath *path) { if (path == (plPath *)NULL) return; if (path->type == PATH_SEGMENT_LIST && path->segments_len > 0) /* number of slots allocated */ free (path->segments); path->segments = (plPathSegment *)NULL; path->segments_len = 0; path->type = PATH_SEGMENT_LIST; /* restore to default */ path->num_segments = 0; path->primitive = false; path->llx = DBL_MAX; path->lly = DBL_MAX; path->urx = -(DBL_MAX); path->ury = -(DBL_MAX); } void _add_moveto (plPath *path, plPoint p) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; /* empty, so allocate a segment buffer */ path->segments = (plPathSegment *) _pl_xmalloc (DATAPOINTS_BUFSIZ * sizeof(plPathSegment)); path->segments_len = DATAPOINTS_BUFSIZ; path->segments[0].type = S_MOVETO; path->segments[0].p = p; path->num_segments = 1; path->llx = p.x; path->lly = p.y; path->urx = p.x; path->ury = p.y; } void _add_line (plPath *path, plPoint p) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; if (path->num_segments == 0) /* empty, so allocate a segment buffer */ { path->segments = (plPathSegment *) _pl_xmalloc (DATAPOINTS_BUFSIZ * sizeof(plPathSegment)); path->segments_len = DATAPOINTS_BUFSIZ; } if (path->num_segments == path->segments_len) /* full, so reallocate */ { path->segments = (plPathSegment *) _pl_xrealloc (path->segments, 2 * path->segments_len * sizeof(plPathSegment)); path->segments_len *= 2; } path->segments[path->num_segments].type = S_LINE; path->segments[path->num_segments].p = p; path->num_segments++; path->llx = DMIN(path->llx, p.x); path->lly = DMIN(path->lly, p.y); path->urx = DMAX(path->urx, p.x); path->ury = DMAX(path->ury, p.y); } void _add_closepath (plPath *path) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; if (path->num_segments == 0) /* meaningless */ return; if (path->num_segments == path->segments_len) /* full, so reallocate */ { path->segments = (plPathSegment *) _pl_xrealloc (path->segments, 2 * path->segments_len * sizeof(plPathSegment)); path->segments_len *= 2; } path->segments[path->num_segments].type = S_CLOSEPATH; path->segments[path->num_segments].p = path->segments[0].p; path->num_segments++; } void _add_bezier2 (plPath *path, plPoint pc, plPoint p) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; if (path->num_segments == 0) /* empty, so allocate a segment buffer */ { path->segments = (plPathSegment *) _pl_xmalloc (DATAPOINTS_BUFSIZ * sizeof(plPathSegment)); path->segments_len = DATAPOINTS_BUFSIZ; } if (path->num_segments == path->segments_len) /* full, so reallocate */ { path->segments = (plPathSegment *) _pl_xrealloc (path->segments, 2 * path->segments_len * sizeof(plPathSegment)); path->segments_len *= 2; } path->segments[path->num_segments].type = S_QUAD; path->segments[path->num_segments].p = p; path->segments[path->num_segments].pc = pc; path->num_segments++; } void _add_bezier3 (plPath *path, plPoint pc, plPoint pd, plPoint p) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; if (path->num_segments == 0) /* empty, so allocate a segment buffer */ { path->segments = (plPathSegment *) _pl_xmalloc (DATAPOINTS_BUFSIZ * sizeof(plPathSegment)); path->segments_len = DATAPOINTS_BUFSIZ; } if (path->num_segments == path->segments_len) /* full, so reallocate */ { path->segments = (plPathSegment *) _pl_xrealloc (path->segments, 2 * path->segments_len * sizeof(plPathSegment)); path->segments_len *= 2; } path->segments[path->num_segments].type = S_CUBIC; path->segments[path->num_segments].p = p; path->segments[path->num_segments].pc = pc; path->segments[path->num_segments].pd = pd; path->num_segments++; } void _add_arc (plPath *path, plPoint pc, plPoint p1) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; if (path->num_segments == 0) /* empty, so allocate a segment buffer */ { path->segments = (plPathSegment *) _pl_xmalloc (DATAPOINTS_BUFSIZ * sizeof(plPathSegment)); path->segments_len = DATAPOINTS_BUFSIZ; } if (path->num_segments == path->segments_len) /* full, so reallocate */ { path->segments = (plPathSegment *) _pl_xrealloc (path->segments, 2 * path->segments_len * sizeof(plPathSegment)); path->segments_len *= 2; } path->segments[path->num_segments].type = S_ARC; path->segments[path->num_segments].p = p1; path->segments[path->num_segments].pc = pc; path->num_segments++; } void _add_ellarc (plPath *path, plPoint pc, plPoint p1) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; if (path->num_segments == 0) /* empty, so allocate a segment buffer */ { path->segments = (plPathSegment *) _pl_xmalloc (DATAPOINTS_BUFSIZ * sizeof(plPathSegment)); path->segments_len = DATAPOINTS_BUFSIZ; } if (path->num_segments == path->segments_len) /* full, so reallocate */ { path->segments = (plPathSegment *) _pl_xrealloc (path->segments, 2 * path->segments_len * sizeof(plPathSegment)); path->segments_len *= 2; } path->segments[path->num_segments].type = S_ELLARC; path->segments[path->num_segments].p = p1; path->segments[path->num_segments].pc = pc; path->num_segments++; } void _add_box (plPath *path, plPoint p0, plPoint p1, bool clockwise) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; path->type = PATH_BOX; path->p0 = p0; path->p1 = p1; path->clockwise = clockwise; path->llx = DMIN(path->llx, p0.x); path->lly = DMIN(path->lly, p0.y); path->urx = DMAX(path->urx, p0.x); path->ury = DMAX(path->ury, p0.y); path->llx = DMIN(path->llx, p1.x); path->lly = DMIN(path->lly, p1.y); path->urx = DMAX(path->urx, p1.x); path->ury = DMAX(path->ury, p1.y); } void _add_circle (plPath *path, plPoint pc, double radius, bool clockwise) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; path->type = PATH_CIRCLE; path->pc = pc; path->radius = radius; path->clockwise = clockwise; } void _add_ellipse (plPath *path, plPoint pc, double rx, double ry, double angle, bool clockwise) { if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; path->type = PATH_ELLIPSE; path->pc = pc; path->rx = rx; path->ry = ry; path->angle = angle; path->clockwise = clockwise; } /* Draw a polygonal approximation to the circular arc from p0 to p1, with center pc, by calling _fakearc(), which in turn repeatedly calls _add_line(). It is assumed that p0 and p1 are distinct. It is assumed that pc is on the perpendicular bisector of the line segment joining them, and that the graphics cursor is initially located at p0. */ void _add_arc_as_lines (plPath *path, plPoint pc, plPoint p1) { /* starting point */ plPoint p0; /* bisection point of arc, and midpoint of chord */ plPoint pb, pm; /* rotation matrix */ double m[4]; /* other variables */ plVector v, v0, v1; double radius, sagitta; double cross, orientation; /* handcrafted relative chordal deviation table, for this arc */ double custom_chord_table[TABULATED_ARC_SUBDIVISIONS]; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; /* determine starting point */ p0 = path->segments[path->num_segments - 1].p; if (p0.x == p1.x && p0.y == p1.y) /* zero-length arc, draw as zero-length line segment */ _add_line (path, p0); else /* genuine polygonal approximation */ { /* vectors from pc to p0, and pc to p1 */ v0.x = p0.x - pc.x; v0.y = p0.y - pc.y; v1.x = p1.x - pc.x; v1.y = p1.y - pc.y; /* cross product, zero if points are collinear */ cross = v0.x * v1.y - v1.x * v0.y; /* Compute orientation. Note libplot convention: if p0, p1, pc are collinear then arc goes counterclockwise from p0 to p1. */ orientation = (cross >= 0.0 ? 1.0 : -1.0); radius = DIST(pc, p0); /* radius is distance to p0 or p1 */ pm.x = 0.5 * (p0.x + p1.x); /* midpoint of chord from p0 to p1 */ pm.y = 0.5 * (p0.y + p1.y); v.x = p1.x - p0.x; /* chord vector from p0 to p1 */ v.y = p1.y - p0.y; _vscale(&v, radius); pb.x = pc.x + orientation * v.y; /* bisection point of arc */ pb.y = pc.y - orientation * v.x; sagitta = DIST(pb, pm) / radius; /* fill in entries of chordal deviation table for this user-defined sagitta */ _prepare_chord_table (sagitta, custom_chord_table); /* call _fakearc(), using for `rotation' matrix m[] a clockwise or counterclockwise rotation by 90 degrees, depending on orientation */ m[0] = 0.0, m[1] = orientation, m[2] = -orientation, m[3] = 0.0; _fakearc (path, p0, p1, USER_DEFINED_ARC, custom_chord_table, m); } } /* Draw a polygonal approximation to a quarter-ellipse from p0 to p1, by calling _fakearc(), which in turn repeatedly calls _add_line(). pc is the center of the arc, and p0, p1, pc are assumed not to be collinear. It is assumed that the graphics cursor is located at p0 when this function is called. The control triangle for the elliptic arc will have vertices p0, p1, and K = p0 + (p1 - pc) = p1 + (p0 - pc). The arc will pass through p0 and p1, and will be tangent at p0 to the edge from p0 to K, and at p1 to the edge from p1 to K. */ void _add_ellarc_as_lines (plPath *path, plPoint pc, plPoint p1) { plPoint p0; plVector v0, v1; double cross; double m[4]; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; /* determine starting point */ p0 = path->segments[path->num_segments - 1].p; /* vectors from pc to p0, and pc to p1 */ v0.x = p0.x - pc.x; v0.y = p0.y - pc.y; v1.x = p1.x - pc.x; v1.y = p1.y - pc.y; /* cross product */ cross = v0.x * v1.y - v1.x * v0.y; if (FROUND(cross) == 0.0) /* collinear points, draw line segment from p0 to p1 (not quite right, could be bettered) */ _add_line (path, p1); else { /* `rotation' matrix (it maps v0 -> -v1 and v1 -> v0) */ m[0] = - (v0.x * v0.y + v1.x * v1.y) / cross; m[1] = (v0.x * v0.x + v1.x * v1.x) / cross; m[2] = - (v0.y * v0.y + v1.y * v1.y) / cross; m[3] = (v0.x * v0.y + v1.x * v1.y) / cross; /* draw polyline inscribed in the quarter-ellipse */ _fakearc (path, p0, p1, QUARTER_ARC, (double *)NULL, m); } } /* A function that approximates a circular arc by a cubic Bezier. The approximation used is a standard one. E.g., a quarter circle extending from (1,0) to (0,1), with center (0,0), would be approximated by a cubic Bezier with control points (1,KAPPA) and (KAPPA,1). Here KAPPA = (4/3)[sqrt(2)-1] = 0.552284749825, approximately. The cubic Bezier will touch the quarter-circle along the line x=y. For a quarter-circle, the maximum relative error in r as a function of theta is about 2.7e-4. The error in r has the same sign, for all theta. */ /* According to Berthold K. P. Horn , the general formula for KAPPA, for a radius-1 circular arc (not necessary a quarter-circle), KAPPA = (4/3)sqrt[(1-cos H)/(1+cos H)] = (4/3)[1-cos H]/[sin H] = (4/3)[sin H]/[1+cosH] where H is half the angle subtended by the arc. H=45 degrees for a quarter circle. This is the formula we use. */ /* Louis Vosloo points out that for a quarter-circle, the value 0.55228... for KAPPA is, for some purposes, sub-optimal. By dropping the requirement that the quarter-circle and the Bezier touch each other along the symmetry line x=y, one can slightly decrease the maximum relative error. He says 0.5541... is the best possible choice. He doesn't have an improved value of KAPPA for a general arc, though. */ void _add_arc_as_bezier3 (plPath *path, plPoint pc, plPoint p1) { plPoint p0; plVector v0, v1; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; /* determine starting point */ p0 = path->segments[path->num_segments - 1].p; /* vectors to starting, ending points */ v0.x = p0.x - pc.x; v0.y = p0.y - pc.y; v1.x = p1.x - pc.x; v1.y = p1.y - pc.y; if ((v0.x == 0.0 && v0.y == 0.0) || (v1.x == 0.0 && v1.y == 0.0) || (v0.x == v1.x && v0.y == v1.y)) /* degenerate case */ _add_line (path, p1); else /* normal case */ { double oldangle, newangle, anglerange; double cross; int orientation; /* cross product, zero if points are collinear */ cross = v0.x * v1.y - v1.x * v0.y; /* Compute orientation. Note libplot convention: if p0, p1, pc are collinear then arc goes counterclockwise from p0 to p1. */ orientation = (cross >= 0.0 ? 1 : -1); /* compute signed subtended angle */ oldangle = _xatan2 (v0.y, v0.x); newangle = _xatan2 (v1.y, v1.x); anglerange = newangle - oldangle; if (anglerange > M_PI) anglerange -= (2 * M_PI); if (anglerange <= -(M_PI)) anglerange += (2 * M_PI); if (FABS(anglerange) > 0.51 * M_PI) /* subtended angle > 0.51 * pi, so split arc in two and recurse, since Bezier approximation isn't very good for angles much greater than 90 degrees */ { double radius; plPoint pb; plVector v; radius = DIST(pc, p0); /* radius is distance to p0 or p1 */ v.x = p1.x - p0.x; /* chord vector from p0 to p1 */ v.y = p1.y - p0.y; _vscale(&v, radius); pb.x = pc.x + orientation * v.y; /* bisection point of arc */ pb.y = pc.y - orientation * v.x; _add_arc_as_bezier3 (path, pc, pb); _add_arc_as_bezier3 (path, pc, p1); } else /* subtended angle <= 0.51 * pi, so a single Bezier suffices */ { double halfangle, sinhalf, coshalf, kappa; plPoint pc_bezier3, pd_bezier3; halfangle = 0.5 * FABS(anglerange); sinhalf = sin (halfangle); coshalf = cos (halfangle); /* compute kappa using either of two formulae, depending on numerical stability */ if (FABS(sinhalf) < 0.5) kappa = (4.0/3.0) * sinhalf / (1.0 + coshalf); else kappa = (4.0/3.0) * (1.0 - coshalf) / sinhalf; pc_bezier3.x = p0.x - kappa * orientation * v0.y; pc_bezier3.y = p0.y + kappa * orientation * v0.x; pd_bezier3.x = p1.x + kappa * orientation * v1.y; pd_bezier3.y = p1.y - kappa * orientation * v1.x; _add_bezier3 (path, pc_bezier3, pd_bezier3, p1); } } } #define KAPPA_FOR_QUARTER_CIRCLE 0.552284749825 void _add_ellarc_as_bezier3 (plPath *path, plPoint pc, plPoint p1) { plPoint p0, pc_bezier3, pd_bezier3; plVector v0, v1; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; /* determine starting point */ p0 = path->segments[path->num_segments - 1].p; /* vectors to starting, ending points */ v0.x = p0.x - pc.x; v0.y = p0.y - pc.y; v1.x = p1.x - pc.x; v1.y = p1.y - pc.y; /* replace by cubic Bezier, with computed control points */ pc_bezier3.x = p0.x + KAPPA_FOR_QUARTER_CIRCLE * v1.x; pc_bezier3.y = p0.y + KAPPA_FOR_QUARTER_CIRCLE * v1.y; pd_bezier3.x = p1.x + KAPPA_FOR_QUARTER_CIRCLE * v0.x; pd_bezier3.y = p1.y + KAPPA_FOR_QUARTER_CIRCLE * v0.y; _add_bezier3 (path, pc_bezier3, pd_bezier3, p1); } /* Approximate a quadratic Bezier by a polyline: standard deCasteljau bisection algorithm. However, we stop subdividing when an appropriate metric of the quadratic Bezier to be drawn is sufficiently small. If (p0,p1,p2) defines the quadratic Bezier, we require that the length of p0-2*p1+p2 be less than REL_QUAD_FLATNESS times the distance between the endpoints of the original Bezier. */ void _add_bezier2_as_lines (plPath *path, plPoint pc, plPoint p) { plPoint r0[MAX_NUM_BEZIER2_SUBDIVISIONS + 1], r1[MAX_NUM_BEZIER2_SUBDIVISIONS + 1], r2[MAX_NUM_BEZIER2_SUBDIVISIONS + 1]; int level[MAX_NUM_BEZIER2_SUBDIVISIONS + 1]; int n = 0; /* index of top of stack, < MAX_NUM_BEZIER2_SUBDIVISIONS */ int segments_drawn = 0; plPoint p0; double sqdist, max_squared_length; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; /* determine starting point */ p0 = path->segments[path->num_segments - 1].p; /* squared distance between p0 and p */ sqdist = (p.x - p0.x) * (p.x - p0.x) + (p.y - p0.y) * (p.y - p0.y); max_squared_length = REL_QUAD_FLATNESS * REL_QUAD_FLATNESS * sqdist; r0[0] = p0; r1[0] = pc; r2[0] = p; level[0] = 0; while (n >= 0) /* i.e. while stack is nonempty */ { int current_level; plPoint q0, q1, q2; current_level = level[n]; q0 = r0[n]; q1 = r1[n]; q2 = r2[n]; if (current_level >= MAX_NUM_BEZIER2_SUBDIVISIONS) /* to avoid stack overflow, draw as line segment */ { _add_line (path, q2); segments_drawn++; n--; } else /* maybe bisect the Bezier */ { plPoint qq0, qq1; plPoint qqq0; plVector vec1; vec1.x = q0.x - 2 * q1.x + q2.x; vec1.y = q0.y - 2 * q1.y + q2.y; if (vec1.x * vec1.x + vec1.y * vec1.y < max_squared_length) /* very flat Bezier, so draw as line segment */ { _add_line (path, q2); segments_drawn++; n--; } else /* split Bezier into pair and recurse */ /* level[n] >= n is an invariant */ { qq0.x = MIDWAY(q0.x, q1.x); qq0.y = MIDWAY(q0.y, q1.y); qq1.x = MIDWAY(q1.x, q2.x); qq1.y = MIDWAY(q1.y, q2.y); qqq0.x = MIDWAY(qq0.x, qq1.x); qqq0.y = MIDWAY(qq0.y, qq1.y); /* first half, deal with next */ r0[n+1] = q0; r1[n+1] = qq0; r2[n+1] = qqq0; level[n+1] = current_level + 1; /* second half, deal with later */ r0[n] = qqq0; r1[n] = qq1; r2[n] = q2; level[n] = current_level + 1; n++; } } } } /* Approximate a cubic Bezier by a polyline: standard deCasteljau bisection algorithm. However, we stop subdividing when an appropriate metric of the cubic Bezier to be drawn is sufficiently small. If (p0,p1,p2,p3) defines the cubic Bezier, we require that the lengths of p0-2*p1+p2 and p1-2*p2+p3 be less than REL_CUBIC_FLATNESS times the distance between the endpoints of the original Bezier. */ void _add_bezier3_as_lines (plPath *path, plPoint pc, plPoint pd, plPoint p) { plPoint r0[MAX_NUM_BEZIER3_SUBDIVISIONS + 1], r1[MAX_NUM_BEZIER3_SUBDIVISIONS + 1], r2[MAX_NUM_BEZIER3_SUBDIVISIONS + 1], r3[MAX_NUM_BEZIER3_SUBDIVISIONS + 1]; int level[MAX_NUM_BEZIER3_SUBDIVISIONS + 1]; int n = 0; /* index of top of stack, < MAX_NUM_BEZIER3_SUBDIVISIONS */ int segments_drawn = 0; plPoint p0; double sqdist, max_squared_length; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments == 0) return; /* determine starting point */ p0 = path->segments[path->num_segments - 1].p; /* squared distance between p0 and p */ sqdist = (p.x - p0.x) * (p.x - p0.x) + (p.y - p0.y) * (p.y - p0.y); max_squared_length = REL_CUBIC_FLATNESS * REL_CUBIC_FLATNESS * sqdist; r0[0] = p0; r1[0] = pc; r2[0] = pd; r3[0] = p; level[0] = 0; while (n >= 0) /* i.e. while stack is nonempty */ { int current_level; plPoint q0, q1, q2, q3; current_level = level[n]; q0 = r0[n]; q1 = r1[n]; q2 = r2[n]; q3 = r3[n]; if (current_level >= MAX_NUM_BEZIER3_SUBDIVISIONS) /* draw line segment, to avoid stack overflow */ { _add_line (path, q3); segments_drawn++; n--; } else /* maybe bisect the Bezier */ { plPoint qq0, qq1, qq2; plPoint qqq0, qqq1; plPoint qqqq0; plVector vec1, vec2; vec1.x = q0.x - 2 * q1.x + q2.x; vec1.y = q0.y - 2 * q1.y + q2.y; vec2.x = q1.x - 2 * q2.x + q3.x; vec2.y = q1.y - 2 * q2.y + q3.y; if (vec1.x * vec1.x + vec1.y * vec1.y < max_squared_length && vec2.x * vec2.x + vec2.y * vec2.y < max_squared_length) /* very flat Bezier, so draw as line segment */ { _add_line (path, q3); segments_drawn++; n--; } else /* split Bezier into pair and recurse */ /* level[n] >= n is an invariant */ { qq0.x = MIDWAY(q0.x, q1.x); qq0.y = MIDWAY(q0.y, q1.y); qq1.x = MIDWAY(q1.x, q2.x); qq1.y = MIDWAY(q1.y, q2.y); qq2.x = MIDWAY(q2.x, q3.x); qq2.y = MIDWAY(q2.y, q3.y); qqq0.x = MIDWAY(qq0.x, qq1.x); qqq0.y = MIDWAY(qq0.y, qq1.y); qqq1.x = MIDWAY(qq1.x, qq2.x); qqq1.y = MIDWAY(qq1.y, qq2.y); qqqq0.x = MIDWAY(qqq0.x, qqq1.x); qqqq0.y = MIDWAY(qqq0.y, qqq1.y); /* first half, deal with next */ level[n+1] = current_level + 1; r0[n+1] = q0; r1[n+1] = qq0; r2[n+1] = qqq0; r3[n+1] = qqqq0; /* second half, deal with later */ level[n] = current_level + 1; r0[n] = qqqq0; r1[n] = qqq1; r2[n] = qq2; r3[n] = q3; n++; } } } } void _add_box_as_lines (plPath *path, plPoint p0, plPoint p1, bool clockwise) { bool x_move_is_first; plPoint newpoint; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; _add_moveto (path, p0); /* if counterclockwise, would first pen motion be in x direction? */ x_move_is_first = ((p1.x >= p0.x && p1.y >= p0.y) || (p1.x < p0.x && p1.y < p0.y) ? true : false); if (clockwise) /* take complement */ x_move_is_first = (x_move_is_first == true ? false : true); if (x_move_is_first) { newpoint.x = p1.x; newpoint.y = p0.y; } else { newpoint.x = p0.x; newpoint.y = p1.y; } _add_line (path, newpoint); _add_line (path, p1); if (x_move_is_first) { newpoint.x = p0.x; newpoint.y = p1.y; } else { newpoint.x = p1.x; newpoint.y = p0.y; } _add_line (path, newpoint); _add_line (path, p0); path->primitive = true; /* flag as flattened primitive */ } void _add_ellipse_as_bezier3s (plPath *path, plPoint pc, double rx, double ry, double angle, bool clockwise) { plPoint startpoint, newpoint; double theta, costheta, sintheta; double xc, yc; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; /* draw ellipse by drawing four elliptic arcs */ theta = (M_PI / 180.0) * angle; /* convert to radians */ costheta = cos (theta); sintheta = sin (theta); xc = pc.x; yc = pc.y; startpoint.x = xc + rx * costheta; startpoint.y = yc + rx * sintheta; _add_moveto (path, startpoint); if (clockwise) { newpoint.x = xc + ry * sintheta; newpoint.y = yc - ry * costheta; } else { newpoint.x = xc - ry * sintheta; newpoint.y = yc + ry * costheta; } _add_ellarc_as_bezier3 (path, pc, newpoint); newpoint.x = xc - rx * costheta; newpoint.y = yc - rx * sintheta; _add_ellarc_as_bezier3 (path, pc, newpoint); if (clockwise) { newpoint.x = xc - ry * sintheta; newpoint.y = yc + ry * costheta; } else { newpoint.x = xc + ry * sintheta; newpoint.y = yc - ry * costheta; } _add_ellarc_as_bezier3 (path, pc, newpoint); _add_ellarc_as_bezier3 (path, pc, startpoint); path->primitive = true; /* flag as flattened primitive */ } void _add_ellipse_as_ellarcs (plPath *path, plPoint pc, double rx, double ry, double angle, bool clockwise) { plPoint startpoint, newpoint; double theta, costheta, sintheta; double xc, yc; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; /* draw ellipse by drawing four elliptic arcs */ theta = (M_PI / 180.0) * angle; /* convert to radians */ costheta = cos (theta); sintheta = sin (theta); xc = pc.x; yc = pc.y; startpoint.x = xc + rx * costheta; startpoint.y = yc + rx * sintheta; _add_moveto (path, startpoint); if (clockwise) { newpoint.x = xc + ry * sintheta; newpoint.y = yc - ry * costheta; } else { newpoint.x = xc - ry * sintheta; newpoint.y = yc + ry * costheta; } _add_ellarc (path, pc, newpoint); newpoint.x = xc - rx * costheta; newpoint.y = yc - rx * sintheta; _add_ellarc (path, pc, newpoint); if (clockwise) { newpoint.x = xc - ry * sintheta; newpoint.y = yc + ry * costheta; } else { newpoint.x = xc + ry * sintheta; newpoint.y = yc - ry * costheta; } _add_ellarc (path, pc, newpoint); _add_ellarc (path, pc, startpoint); path->primitive = true; /* flag as flattened primitive */ } void _add_ellipse_as_lines (plPath *path, plPoint pc, double rx, double ry, double angle, bool clockwise) { plPoint startpoint, newpoint; double theta, costheta, sintheta; double xc, yc; if (path == (plPath *)NULL) return; if (path->type != PATH_SEGMENT_LIST || path->num_segments > 0) return; /* draw ellipse by drawing four fake elliptic arcs */ theta = (M_PI / 180.0) * angle; /* convert to radians */ costheta = cos (theta); sintheta = sin (theta); xc = pc.x; yc = pc.y; startpoint.x = xc + rx * costheta; startpoint.y = yc + rx * sintheta; _add_moveto (path, startpoint); if (clockwise) { newpoint.x = xc + ry * sintheta; newpoint.y = yc - ry * costheta; } else { newpoint.x = xc - ry * sintheta; newpoint.y = yc + ry * costheta; } _add_ellarc_as_lines (path, pc, newpoint); newpoint.x = xc - rx * costheta; newpoint.y = yc - rx * sintheta; _add_ellarc_as_lines (path, pc, newpoint); if (clockwise) { newpoint.x = xc - ry * sintheta; newpoint.y = yc + ry * costheta; } else { newpoint.x = xc + ry * sintheta; newpoint.y = yc - ry * costheta; } _add_ellarc_as_lines (path, pc, newpoint); _add_ellarc_as_lines (path, pc, startpoint); path->primitive = true; /* flag as flattened primitive */ } void _add_circle_as_bezier3s (plPath *path, plPoint pc, double radius, bool clockwise) { if (path == (plPath *)NULL) return; _add_ellipse_as_bezier3s (path, pc, radius, radius, 0.0, clockwise); path->primitive = true; /* flag as flattened primitive */ } void _add_circle_as_ellarcs (plPath *path, plPoint pc, double radius, bool clockwise) { if (path == (plPath *)NULL) return; _add_ellipse_as_ellarcs (path, pc, radius, radius, 0.0, clockwise); path->primitive = true; /* flag as flattened primitive */ } void _add_circle_as_lines (plPath *path, plPoint pc, double radius, bool clockwise) { if (path == (plPath *)NULL) return; _add_ellipse_as_lines (path, pc, radius, radius, 0.0, clockwise); path->primitive = true; /* flag as flattened primitive */ } /* The _fakearc() subroutine below, which is called by _add_arc_as_lines() and _add_ellarc_as_lines(), contains a remote descendent of the arc-drawing algorithm of Ken Turkowski described in Graphics Gems V. His algorithm is a recursive circle subdivision algorithm, which relies on the fact that if s and s' are the (chordal deviation)/radius associated to (respectively) an arc and a half-arc, then s' is approximately equal to s/4. The exact formula is s' = 1 - sqrt (1 - s/2), which applies for all s in the meaningful range, i.e. 0 <= s <= 2. Ken's algorithm rotates the chord of an arc by 90 degrees and scales it to have length s'. The resulting vector will be the chordal deviation vector of the arc, which gives the midpoint of the arc, and hence the endpoints of each half of the arc. So this drawing technique is recursive. The problem with this approach is that scaling a vector to a specified length requires a square root, so there are two square roots in each subdivision step. One can attempt to remove one of them by noticing that the chord half-length h always satisfies h = sqrt(s * (2-s))). So one can rotate the chord vector by 90 degrees, and multiply its length by s/2h, i.e., s/2sqrt(s * (2-s)), to get the chordal deviation vector. This factor still includes a square root though. Also one still needs to compute a square root in order to proceed from one subdivision step to the next, i.e. to compute s' from s. We can get around the square root problem by drawing only circular arcs with subtended angle of 90 degrees (quarter-circles), or elliptic arcs that are obtained from such quarter-circles by affine transformations (so-called quarter-ellipses). To draw the latter, we need only replace the 90 degree rotation matrix mentioned above by an affine transformation that maps v0->-v1 and v1->v0, where v0 = p0 - pc and v1 = p1 - pc are the displacement vectors from the center of the ellipse to the endpoints of the arc. If we do this, we get an elliptic arc with p0 and p1 as endpoints. The vectors v0 and v1 are said lie along conjugate diameters of the quarter-ellipse. So for drawing quarter-ellipses, the only initial value of s we need to consider is the one for a quarter-circle, which is 1-sqrt(1/2). The successive values of s/2h that will be encountered, after each bisection, are pre-computed and stored in a lookup table, found in g_arc.h. This approach lets us avoid, completely, any computation of square roots during the drawing of quarter-circles and quarter-ellipses. The only square roots we need are precomputed. We don't need any floating point divisions in the main loop either. Of course for other angles than 90 degrees, we precompute a chord table and pass it to _fakearc(). The implementation below does not use recursion (we use a local array, rather than the call stack, to store the sequence of generated points). */ static void _fakearc (plPath *path, plPoint p0, plPoint p1, int arc_type, const double *custom_chord_table, const double m[4]) { plPoint p[NUM_ARC_SUBDIVISIONS + 1], q[NUM_ARC_SUBDIVISIONS + 1]; int level[NUM_ARC_SUBDIVISIONS + 1]; int n = 0; /* index of top of stack, < NUM_ARC_SUBDIVISIONS */ int segments_drawn = 0; const double *our_chord_table; if (arc_type == USER_DEFINED_ARC) our_chord_table = custom_chord_table; else /* custom_chord_table arg ignored */ our_chord_table = _chord_table[arc_type]; p[0] = p0; q[0] = p1; level[0] = 0; while (n >= 0) /* i.e. while stack is nonempty */ { if (level[n] >= NUM_ARC_SUBDIVISIONS) { /* draw line segment */ _add_line (path, q[n]); segments_drawn++; n--; } else /* bisect line segment */ { plVector v; plPoint pm, pb; v.x = q[n].x - p[n].x; /* chord = line segment from p[n] to q[n] */ v.y = q[n].y - p[n].y; pm.x = p[n].x + 0.5 * v.x; /* midpoint of chord */ pm.y = p[n].y + 0.5 * v.y; /* Compute bisection point. If m=[0 1 -1 0] this just rotates the chord clockwise by 90 degrees, and scales it to yield the chordal deviation vector, which is used as an offset. */ pb.x = pm.x + our_chord_table[level[n]] * (m[0] * v.x + m[1] * v.y); pb.y = pm.y + our_chord_table[level[n]] * (m[2] * v.x + m[3] * v.y); /* replace line segment by pair; level[n] >= n is an invariant */ p[n+1] = p[n]; q[n+1] = pb; /* first half, deal with next */ level[n+1] = level[n] + 1; p[n] = pb; q[n] = q[n]; /* second half, deal with later */ level[n] = level[n] + 1; n++; } } } /* prepare_chord_table() computes the list of chordal deviation factors that _fakearc() needs when it is employed to draw a circular arc of subtended angle other than the default angles it supports */ static void _prepare_chord_table (double sagitta, double custom_chord_table[TABULATED_ARC_SUBDIVISIONS]) { double half_chord_length; int i; half_chord_length = sqrt ( sagitta * (2.0 - sagitta) ); for (i = 0; i < TABULATED_ARC_SUBDIVISIONS; i++) { custom_chord_table[i] = 0.5 * sagitta / half_chord_length; sagitta = 1.0 - sqrt (1.0 - 0.5 * sagitta); half_chord_length = 0.5 * half_chord_length / (1.0 - sagitta); } } /* Flatten a simple path into a path of segment list type, consisting only of line segments. As supplied, the path may be perfectly general: it may be a segment list (not all segments necessarily being line segments), or be a closed primitive (box/circle/ellipse). If supplied path already consists only of line segments (with an initial moveto and possibly a final closepath), it is returned unchanged; this can be tested for by comparing pointers for equality. If a new path is returned, it must be freed with _delete_plPath(). */ plPath * _flatten_path (const plPath *path) { plPath *newpath; if (path == (plPath *)NULL) return (plPath *)NULL; switch (path->type) { case PATH_SEGMENT_LIST: { bool do_flatten = false; int i; for (i = 0; i < path->num_segments; i++) { if (path->segments[i].type != S_MOVETO && path->segments[i].type != S_LINE && path->segments[i].type != S_CLOSEPATH) { do_flatten = true; break; } } if (do_flatten == false) newpath = (plPath *)path; /* just return original path */ else { newpath = _new_plPath (); for (i = 0; i < path->num_segments; i++) { switch ((int)(path->segments[i].type)) { case (int)S_MOVETO: _add_moveto (newpath, path->segments[i].p); break; case (int)S_LINE: _add_line (newpath, path->segments[i].p); break; case (int)S_CLOSEPATH: _add_closepath (newpath); break; /* now, the types of segment we flatten: */ case (int)S_ARC: _add_arc_as_lines (newpath, path->segments[i].pc, path->segments[i].p); break; case (int)S_ELLARC: _add_ellarc_as_lines (newpath, path->segments[i].pc, path->segments[i].p); break; case (int)S_QUAD: _add_bezier2_as_lines (newpath, path->segments[i].pc, path->segments[i].p); break; case (int)S_CUBIC: _add_bezier3_as_lines (newpath, path->segments[i].pc, path->segments[i].pd, path->segments[i].p); break; default: /* shouldn't happen */ break; } } } break; } case PATH_CIRCLE: newpath = _new_plPath (); _add_circle_as_lines (newpath, path->pc, path->radius, path->clockwise); break; case PATH_ELLIPSE: newpath = _new_plPath (); _add_ellipse_as_lines (newpath, path->pc, path->rx, path->ry, path->angle, path->clockwise); break; case PATH_BOX: newpath = _new_plPath (); _add_box_as_lines (newpath, path->p0, path->p1, path->clockwise); break; default: /* shouldn't happen */ newpath = _new_plPath (); break; } return newpath; } /**********************************************************************/ /* The code below exports the _merge_paths() function, which munges an array of plPath objects and returns a new array. Heuristic "path-merging" of this sort is performed in g_endpath.c when filling a compound path (i.e., a multi-plPath path), if the output driver supports only the filling of single plPaths. That is the case for nearly all of libplot's output drivers. `Child' paths are merged into their `parents', so each location in the returned array where a child was present will contain NULL. Also, any non-child that had no children will be returned without modification. You should take this into account when freeing the returned array of plPaths. Only the elements that are (1) non-NULL, and (2) differ from the corresponding elements in the originally passed array should have _delete_plPath() invoked on them. The array as a whole may be deallocated by calling free(). The _merge_paths() function was inspired by a similar function in Wolfgang Glunz's pstoedit, which was originally written by Burkhard Plaum . */ /* Note: a well-formed plPath has the form: { moveto { line | arc }* { closepath }? } The _merge_paths function was written to merge _closed_ plPath's, i.e. ones whose endpoint is the same as the startpoint (possibly implicitly, i.e. closepath is allowed). However, libplot applies it to open paths too, in which case an `implicit closepath' is added to the path to close it. NOTE: The current release of libplot does not yet support `closepath' segments at a higher level. So we regard a pass-in plPath as `closed' if its last defining vertex is the same as the first. THIS CONVENTION WILL GO AWAY. */ /* ad hoc structure for an annotated plPath, in particular one that has been flattened into line segments and annotated; used only in this file, for merging purposes */ typedef struct subpath_struct { plPathSegment *segments; /* segment array */ int num_segments; /* number of segments */ struct subpath_struct **parents; /* array of pointers to possible parents */ struct subpath_struct *parent; /* pointer to parent path */ struct subpath_struct **children; /* array of pointers to child paths */ int num_children; /* number of children */ int num_outside; /* number of subpaths outside this one */ double llx, lly, urx, ury; /* bounding box of the subpath */ bool inserted; /* subpath has been inserted into result? */ } subpath; /* forward references */ /* 0. ctors, dtors */ static subpath * new_subpath (void); static subpath ** new_subpath_array (int n); static void delete_subpath (subpath *s); static void delete_subpath_array (subpath **s, int n); /* 1. functions that act on a subpath, i.e. an `annotated path' */ static bool is_inside_of (const subpath *s, const subpath *other); static double _cheap_lower_bound_on_distance (const subpath *path1, const subpath *path2); static void linearize_subpath (subpath *s); static void read_into_subpath (subpath *s, const plPath *path); /* 2. miscellaneous */ static void find_parents_in_subpath_list (subpath **annotated_paths, int num_paths); static void insert_subpath (plPathSegment *parent_segments, const plPathSegment *child_segments, int parent_size, int child_size, int parent_index, int child_index); static void _compute_closest (const plPathSegment *p1, const plPathSegment *p2, int size1, int size2, double *distance, int *index1, int *index2); /**********************************************************************/ /* ctor for subpath class */ static subpath * new_subpath (void) { subpath *s; s = (subpath *)_pl_xmalloc (sizeof (subpath)); s->segments = (plPathSegment *)NULL; s->num_segments = 0; s->parents = (subpath **)NULL; s->parent = (subpath *)NULL; s->children = (subpath **)NULL; s->num_children = 0; s->num_outside = 0; s->llx = DBL_MAX; s->lly = DBL_MAX; s->urx = -DBL_MAX; s->ury = -DBL_MAX; s->inserted = false; return s; } /* corresponding ctor for a subpath array */ static subpath ** new_subpath_array (int n) { int i; subpath **s; s = (subpath **)_pl_xmalloc (n * sizeof (subpath *)); for (i = 0; i < n; i++) s[i] = new_subpath (); return s; } /* dtor for subpath class */ static void delete_subpath (subpath *s) { if (s) { if (s->segments) free (s->segments); if (s->children) free (s->children); if (s->parents) free (s->parents); free (s); } } /* corresponding dtor for a subpath array */ static void delete_subpath_array (subpath **s, int n) { int i; if (s) { for (i = 0; i < n; i++) delete_subpath (s[i]); free (s); } } /* replace every segment in a subpath by a lineto (invoked only on a child subpath, i.e. a subpath with an identified parent) */ static void linearize_subpath (subpath *s) { /* replace first segment (moveto) with a lineto */ s->segments[0].type = S_LINE; /* if final segment is a closepath, also replace with a lineto, back to point #0 */ if (s->segments[s->num_segments-1].type == S_CLOSEPATH) { s->segments[s->num_segments-1].type = S_LINE; s->segments[s->num_segments-1].p = s->segments[0].p; } } /* Read a sequence of plPathSegments from a plPath into a previously empty annotated subpath. (This is called only after the plPath has been flattened, so the segments include no arcs.) Because the way in which _merge_paths() is currently called in libplot, we need to handle the possibility that the plPath may not be closed. If it isn't, we add a closepath. At present, we allow a final lineto to the start vertex to serve the same purpose. THIS IS A LIBPLOT CONVENTION THAT WILL GO AWAY. */ static void read_into_subpath (subpath *s, const plPath *path) { bool need_to_close = false; int i; /* sanity check */ if (path->type != PATH_SEGMENT_LIST) return; /* allocate space for segment array of subpath; add 1 extra slot for manual closure, if needed */ s->segments = (plPathSegment *)_pl_xmalloc((path->num_segments + 1) * sizeof (plPathSegment)); s->num_segments = path->num_segments; /* sanity check */ if (path->num_segments == 0) return; /* Is this path closed? If not, we'll close manually the annotated path that we'll construct. WE CURRENTLY TREAT FINAL = INITIAL AS INDICATING CLOSURE. */ if (path->segments[path->num_segments - 1].type != S_CLOSEPATH && (path->segments[path->num_segments - 1].p.x != path->segments[0].p.x || path->segments[path->num_segments - 1].p.y != path->segments[0].p.y)) need_to_close = true; /* copy the segments, updating bounding box to take each juncture point into account */ for (i = 0; i < path->num_segments; i++) { plPathSegment e; e = path->segments[i]; s->segments[i] = e; if (e.p.x < s->llx) s->llx = e.p.x; if (e.p.y < s->lly) s->lly = e.p.y; if (e.p.x > s->urx) s->urx = e.p.x; if (e.p.y > s->ury) s->ury = e.p.y; } if (need_to_close) { #if 0 s->segments[path->num_segments].type = S_CLOSEPATH; #else /* currently, use line segment instead of closepath */ s->segments[path->num_segments].type = S_LINE; #endif s->segments[path->num_segments].p = path->segments[0].p; s->num_segments++; } } /* check if a subpath is inside another subpath */ static bool is_inside_of (const subpath *s, const subpath *other) { int inside = 0; int outside = 0; int i; /* if bbox fails to lie inside the other's bbox, false */ if (!((s->llx >= other->llx) && (s->lly >= other->lly) && (s->urx <= other->urx) && (s->ury <= other->ury))) return false; /* otherwise, check all juncture points */ for (i = 0; i < s->num_segments; i++) { bool point_is_inside; if (s->segments[i].type == S_CLOSEPATH) /* should have i = num_segments - 1, no associated juncture point */ continue; /* Check if the vertex s->segments[i].p is inside `other'. Could be done in a separate function, but we inline it for speed. */ { /* These two factors should be small positive floating-point numbers. They should preferably be incommensurate, to minimize the probability of a degenerate case occurring: two line segments intersecting at the endpoint of one or the other. */ #define SMALL_X_FACTOR (M_SQRT2 * M_PI) #define SMALL_Y_FACTOR (M_SQRT2 + M_PI) /* argument of the now-inlined function (besides `other') */ plPoint p; /* local variables of the now-inlined function */ int k, crossings; /* (x1,y1) is effectively the point at infinity */ double x1, y1; /* (x2,y2) is specified point */ double x2, y2; /* argument of the now-inlined function (besides `other') */ p = s->segments[i].p; /* (x1,y1) is effectively the point at infinity */ x1 = (DMAX(p.x, other->urx) + SMALL_X_FACTOR * (DMAX(p.x, other->urx) - DMIN(p.x, other->llx))); y1 = (DMAX(p.y, other->ury) + SMALL_Y_FACTOR * (DMAX(p.y, other->ury) - DMIN(p.y, other->lly))); /* (x2,y2) is specified point */ x2 = p.x; y2 = p.y; crossings = 0; for (k = 0; k < other->num_segments; k++) { int j; double x3, y3, x4, y4, det, det1, det2; if (other->segments[k].type == S_CLOSEPATH) /* k > 0 */ { x3 = other->segments[k-1].p.x; y3 = other->segments[k-1].p.y; } else { x3 = other->segments[k].p.x; y3 = other->segments[k].p.y; } j = (k == other->num_segments - 1 ? 0 : k + 1); if (other->segments[j].type == S_CLOSEPATH) continue; x4 = other->segments[j].p.x; y4 = other->segments[j].p.y; /* (x3,y3)-(x4,y4) is a line segment in the closed path */ /* Check whether the line segments (x1,y1)-(x2,y2) and (x3-y3)-(x4,y4) cross each other. System to solve is: [p1 + (p2 - p1) * t1] - [p3 + (p4 - p3) * t2] = 0 i.e. (x2 - x1) * t1 - (x4 - x3) * t2 = x3 - x1; (y2 - y1) * t1 - (y4 - y3) * t2 = y3 - y1; Solutions are: t1 = det1/det t2 = det2/det The line segments cross each other (in their interiors) if 0.0 < t1 < 1.0 and 0.0 < t2 < 1.0 */ det = (x2 - x1) * (-(y4 - y3)) - (-(x4 - x3)) * (y2 - y1); if (det == 0.0) /* line segments are parallel; ignore the degenerate case that they might overlap */ continue; det1 = (x3 - x1) * (-(y4 - y3)) - (-(x4 - x3)) * (y3 - y1); det2 = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1); if ((det<0.0 && (det1>0.0 || det2>0.0 || det10.0 && (det1<0.0 || det2<0.0 || det1>det || det2>det))) /* solution for at least one of t1 and t2 is outside the interval [0,1], so line segments do not cross */ continue; /* We ignore the possibility that t1, t2 are both in the interval [0,1], but (t1 == 0.0) || (t1 == 1.0) || (t2 == 0.0) || (t2 == 1.0). t1 == 0.0 should never happen, if p1 is effectively the point at infinity. So this degenerate case occurs only if the line segment (x1,y1)-(x2,y2) goes through either (x3,y3) or (x4,y4), or the specified point (x2,y2) lies on the line segment (x3,y3)-(x4,y4) that is part of the path. */ crossings++; } /* our determination of whether the point is inside the path; before we inlined this function, this was the return value */ point_is_inside = (crossings & 1) ? true : false; } /* increment inside,outside depending on whether or not the juncture point was inside the other path */ if (point_is_inside) inside++; else outside++; } /* make a democratic decision as to whether the path as a whole is inside the other path */ return (inside > outside ? true : false); } /* Find parent (if any) of each subpath in a list of subpaths. When this is invoked, each subpath should consist of an initial moveto, at least one lineto, and a closepath (not currently enforced). */ static void find_parents_in_subpath_list (subpath **annotated_paths, int num_paths) { int i, j; subpath *parent; /* determine for each subpath the subpaths that are nominally outside it */ for (i = 0; i < num_paths; i++) { annotated_paths[i]->parents = new_subpath_array (num_paths); for (j = 0; j < num_paths; j++) { if (j != i) { if (is_inside_of (annotated_paths[i], annotated_paths[j])) { annotated_paths[i]->parents[annotated_paths[i]->num_outside] = annotated_paths[j]; annotated_paths[i]->num_outside++; } } } } /* Now find the real parent subpaths, i.e. the root subpaths. A subpath is a parent subpath if the number of nominally-outside subpaths is even, and is a child subpath only if the number is odd. An odd number, together with a failure to find a suitable potential parent, will flag a path as an isolate: technically a parent, but without children. */ for (i = 0; i < num_paths; i++) { if ((annotated_paths[i]->num_outside & 1) == 0) /* an even number of subpaths outside, definitely a parent (i.e. doesn't have a parent itself, may or may not have children) */ { /* allocate space for children (if any) */ annotated_paths[i]->children = new_subpath_array (num_paths); } } /* now determine which are children, and link them to their parents */ for (i = 0; i < num_paths; i++) { if ((annotated_paths[i]->num_outside & 1) == 0) /* even number outside, definitely a parent subpath (whether it has children remains to be determined) */ continue; else /* odd number outside, possibly a child, so search linearly through possible parents until we get a hit; if so, this is a child; if not, this is an isolate (classed as a parent) */ { for (j = 0; j < annotated_paths[i]->num_outside; j++) { if (annotated_paths[i]->num_outside == annotated_paths[i]->parents[j]->num_outside + 1) /* number outside is one more than the number outside a potential parent; flag as a child, and add it to the parent's child list */ { parent = annotated_paths[i]->parents[j]; annotated_paths[i]->parent = parent; /* give it a parent */ parent->children[parent->num_children] = annotated_paths[i]; parent->num_children++; break; } } } } } /* Compute closest vertices in two paths. Indices of closest vertices, and (squared) distance between them, are returned via pointers. This is invoked in _merge_paths() only on paths that have been flattened, and have had the initial moveto and the optional final closepath replaced by line segments. So when this is called, each segment type is S_LINE. */ static void _compute_closest (const plPathSegment *p1, const plPathSegment *p2, int size1, int size2, double *distance, int *index1, int *index2) { int best_i = 0, best_j = 0; /* keep compiler happy */ double best_distance = DBL_MAX; int ii, jj; for (ii = 0; ii < size1; ii++) { plPoint point1; point1 = p1[ii].p; for (jj = 0; jj < size2; jj++) { double tmp1, tmp2, distance; plPoint point2; point2 = p2[jj].p; tmp1 = point1.x - point2.x; tmp2 = point1.y - point2.y; distance = tmp1 * tmp1 + tmp2 * tmp2; if (distance < best_distance) { best_distance = distance; best_i = ii; best_j = jj; } } } /* return the three quantities */ *distance = best_distance; *index1 = best_i; *index2 = best_j; } /* Compute a cheap lower bound on the (squared) distance between two subpaths, by looking at their bounding boxes. */ static double _cheap_lower_bound_on_distance (const subpath *path1, const subpath *path2) { double xdist = 0.0, ydist = 0.0, dist; if (path1->urx < path2->llx) xdist = path2->llx - path1->urx; else if (path2->urx < path1->llx) xdist = path1->llx - path2->urx; if (path1->ury < path2->lly) ydist = path2->lly - path1->ury; else if (path2->ury < path1->lly) ydist = path1->lly - path2->ury; dist = xdist * xdist + ydist * ydist; return dist; } /* Insert a closed child subpath into a closed parent, by connecting the two, twice, at a specified vertex of the parent path and a specified vertex of the child path. This is the key function invoked by _merge_paths(). Both paths are supplied as segment lists, and all segments are lines; the final segment of each is a line back to the starting point. I.e., for both subpaths, the final vertex duplicates the start vertex. So we ignore the final vertex of the child (but not that of the parent). I.e. if the child vertices are numbered 0..child_size-1, we map the case child_index = child_size-1 to child_index = 0. The new path has parent_size + child_size + 1 vertices, of which the last is a repetition of the first. */ /* INDEX MAP: PARENT -> NEW PATH i --> i (if 0 <= i < parent_index + 1) i --> i+child_size+1 (if parent_index + 1 <= i < parent_size) CHILD -> NEW PATH i --> i+parent_index+child_size-child_index (0 <= i < child_index+1) i --> i+parent_index-child_index+1 (child_index+1 <= i < child_size-1) (0'th element of the child is same as child_size-1 element) NEW PATH CONTENTS 0..parent_index 0..parent_index of PARENT parent_index+1 child_index of CHILD (i.e. ->join) parent_index+2..parent_index+child_size-child_index-1 child_index+1..child_size-2 of CHILD parent_index+child_size-child_index..parent_index+child_size 0..child_index of CHILD parent_index+child_size+1 parent_index of PARENT (i.e. ->join) parent_index+child_size+2..parent_size+child_size parent_index+1..parent_size-1 of PARENT */ /* Macros that map from vertices in the child path and the parent path, to vertices in the merged path. Here the argument `i' is the index in the original path, and each macro evaluates to the index in the merger. */ /* The first macro should not be applied to i=child_size-1; as noted above, that vertex is equivalent to i=0, so apply it to i=0 instead. */ #define CHILD_VERTEX_IN_MERGED_PATH(i,parent_index,parent_size,child_index,child_size) ((i) <= (child_index) ? (i) + (parent_index) + (child_size) - (child_index) : (i) + (parent_index) - (child_index) + 1) #define PARENT_VERTEX_IN_MERGED_PATH(i,parent_index,parent_size,child_index,child_size) ((i) <= (parent_index) ? (i) : (i) + (child_size) + 1) static void insert_subpath (plPathSegment *parent, const plPathSegment *child, int parent_size, int child_size, int parent_index, int child_index) { int i; plPathSegment e1, e2; int src_index; /* map case when joining vertex is final vertex of child to case when it's the 0'th vertex */ if (child_index == child_size - 1) child_index = 0; /* move up: add child_size+1 empty slots to parent path */ for (i = parent_size - 1; i >= parent_index + 1; i--) parent[i + child_size + 1] = parent[i]; /* add a line segment from specified vertex of parent path to specified vertex of child path */ e1 = child[child_index]; e1.type = S_LINE; /* unnecessary */ parent[parent_index + 1] = e1; /* copy vertices of child into parent, looping back to start in child if necessary; note we skip the last (i.e. child_size-1'th) vertex, since the 0'th vertex is the same */ src_index = child_index; for (i = 0; i < child_size - 1; i++) { src_index++; if (src_index == child_size - 1) src_index = 0; parent[parent_index + 2 + i] = child[src_index]; } /* add a line segment back from specified vertex of child path to specified vertex of parent path */ e2 = parent[parent_index]; e2.type = S_LINE; parent[parent_index + child_size + 1] = e2; } /* The key function exported by this module, which is used by libplot for filling compound paths. */ plPath ** _merge_paths (const plPath **paths, int num_paths) { int i; subpath **annotated_paths; plPath **flattened_paths; plPath **merged_paths; /* flatten every path to a list of line segments (some paths may come back unaltered; will be able to compare pointers to check for that) */ flattened_paths = (plPath **)_pl_xmalloc (num_paths * sizeof(plPath *)); for (i = 0; i < num_paths; i++) { flattened_paths[i] = _flatten_path (paths[i]); #ifdef DEBUG fprintf (stderr, "path %d: %d segments, flattened to %d segments\n", i, paths[i]->num_segments, flattened_paths[i]->num_segments); #endif } /* Copy each flattened path into a corresponding annotated path (`subpath'). Manual closure, if necessary (see above) is performed, i.e. we always add a final closepath to close the path. At this stage bounding boxes are computed. */ annotated_paths = new_subpath_array (num_paths); for (i = 0; i < num_paths; i++) read_into_subpath (annotated_paths[i], flattened_paths[i]); /* Flattened paths no longer needed, so delete them carefully (some may be the same as the original paths, due to _flatten_path() having simply returned its argument) */ for (i = 0; i < num_paths; i++) if (flattened_paths[i] != paths[i]) _delete_plPath (flattened_paths[i]); /* determine which subpaths are parents, children */ find_parents_in_subpath_list (annotated_paths, num_paths); /* in each child, replace each moveto/closepath by a lineto */ for (i = 0; i < num_paths; i++) if (annotated_paths[i]->parent != (subpath *)NULL) /* child path */ linearize_subpath (annotated_paths[i]); /* create array of merged paths: parent paths will have child paths merged into them, and child paths won't appear */ /* allocate space for new array, to be returned */ merged_paths = (plPath **)_pl_xmalloc (num_paths * sizeof(plPath *)); for (i = 0; i < num_paths; i++) { int j, k, num_segments_in_merged_path; subpath *parent; plPath *merged_path; double *parent_to_child_distances; int *child_best_indices, *parent_best_indices; if (annotated_paths[i]->parent != (subpath *)NULL) /* child path; original path will be merged into parent */ { merged_paths[i] = (plPath *)NULL; continue; } if (annotated_paths[i]->num_children == 0) /* no parent, but no children either, so no merging done; in output path array, place original unflattened path */ { merged_paths[i] = (plPath *)paths[i]; continue; } /* this path must be a parent, with one or more children to be merged into it; so create new empty `merged path' with segments array that will hold it, and the merged-in children */ parent = annotated_paths[i]; num_segments_in_merged_path = parent->num_segments; for (j = 0; j < parent->num_children; j++) num_segments_in_merged_path += (parent->children[j]->num_segments + 1); merged_path = _new_plPath (); merged_path->segments = (plPathSegment *)_pl_xmalloc(num_segments_in_merged_path * sizeof (plPathSegment)); merged_path->num_segments = 0; merged_path->segments_len = num_segments_in_merged_path; /* copy parent path into new empty path, i.e. initialize the merged path */ for (j = 0; j < parent->num_segments; j++) merged_path->segments[j] = parent->segments[j]; merged_path->num_segments = parent->num_segments; /* Create temporary storage for `closest vertex pairs' and inter-path distances. We first compute the shortest distance between each child path. We also keep track of the shortest distance between each child and the merged path being constructed, and update it when any child is added. */ parent_to_child_distances = (double *)_pl_xmalloc(parent->num_children * sizeof (double)); parent_best_indices = (int *)_pl_xmalloc(parent->num_children * sizeof (int)); child_best_indices = (int *)_pl_xmalloc(parent->num_children * sizeof (int)); /* compute closest vertices between merged path (i.e., right now, the parent) and any child; these arrays will be updated when any child is inserted into the merged path */ for (j = 0; j < parent->num_children; j++) _compute_closest (parent->segments, parent->children[j]->segments, parent->num_segments, parent->children[j]->num_segments, &(parent_to_child_distances[j]), &(parent_best_indices[j]), &(child_best_indices[j])); for (k = 0; k < parent->num_children; k++) /* insert a child (the closest remaining one!) into the built-up merged path; and flag the child as having been inserted so that we don't pay attention to it thereafter */ { double min_distance; int closest = 0; /* keep compiler happy */ double *new_parent_to_child_distances; int *new_child_best_indices, *new_parent_best_indices; /* allocate storage for arrays that will be used to update the three abovementioned arrays, with each pass through the loop */ new_parent_to_child_distances = (double *)_pl_xmalloc(parent->num_children * sizeof (double)); new_parent_best_indices = (int *)_pl_xmalloc(parent->num_children * sizeof (int)); new_child_best_indices = (int *)_pl_xmalloc(parent->num_children * sizeof (int)); /* initially, they're the same as the current arrays */ for (j = 0; j < parent->num_children; j++) { new_parent_to_child_distances[j] = parent_to_child_distances[j]; new_parent_best_indices[j] = parent_best_indices[j]; new_child_best_indices[j] = child_best_indices[j]; } /* find closest child to merged path, which has not yet been inserted */ min_distance = DBL_MAX; for (j = 0; j < parent->num_children; j++) { if (parent->children[j]->inserted) /* ignore this child */ continue; if (parent_to_child_distances[j] < min_distance) { closest = j; min_distance = parent_to_child_distances[j]; } } /* closest remaining child has index `closest'; it will be inserted into the current merged path */ /* loop over all children, skipping inserted ones and also skipping `closest', the next child to be inserted */ for (j = 0; j < parent->num_children; j++) { double inter_child_distance; int inter_child_best_index1, inter_child_best_index2; double lower_bound_on_inter_child_distance; bool compute_carefully; if (parent->children[j]->inserted) /* ignore */ continue; if (j == closest) /* ignore */ continue; /* compute distance (and closest vertex pairs) between `closest' and the j'th child; result is only of interest if the distance is less than parent_to_child_distances[j], so we first compute a cheap lower bound on the result by looking at bounding boxes. */ lower_bound_on_inter_child_distance = _cheap_lower_bound_on_distance (parent->children[j], parent->children[closest]); compute_carefully = (lower_bound_on_inter_child_distance < parent_to_child_distances[j]) ? true : false; if (compute_carefully) /* compute accurate inter-child distance; also which two vertices yield the minimum distance */ _compute_closest (parent->children[j]->segments, parent->children[closest]->segments, parent->children[j]->num_segments, parent->children[closest]->num_segments, &inter_child_distance, &inter_child_best_index1, /* vertex in j */ &inter_child_best_index2); /* in `closest' */ /* fill in j'th element of the new arrays parent_to_child_distances[], parent_best_indices[] and child_best_indices[] so as to take the insertion of the child into account; but we don't update the old arrays until we do the actual insertion */ if (compute_carefully && inter_child_distance < parent_to_child_distances[j]) /* j'th child is nearer to a vertex in `closest', the child to be inserted, than to any vertex in the current merged path, so all three arrays are affected */ { int nearest_index_in_closest_child; new_parent_to_child_distances[j] = inter_child_distance; new_child_best_indices[j] = inter_child_best_index1; nearest_index_in_closest_child = inter_child_best_index2; /* Compute new value of parent_best_indices[j], taking into account that `closest' will be inserted into the merged path, thereby remapping the relevant index in `closest'. The macro doesn't perform correctly if its first arg takes the maximum possible value; so instead, we map that possibility to `0'. See comment above, before the macro definition. */ if (nearest_index_in_closest_child == parent->children[closest]->num_segments - 1) nearest_index_in_closest_child = 0; new_parent_best_indices[j] = CHILD_VERTEX_IN_MERGED_PATH(nearest_index_in_closest_child, parent_best_indices[closest], merged_path->num_segments, child_best_indices[closest], parent->children[closest]->num_segments); } else /* j'th child is nearer to a vertex in the current merged path than to any vertex in `closest', the child to be inserted into the merged path */ { int nearest_index_in_parent; nearest_index_in_parent = parent_best_indices[j]; /* compute new value of parent_best_indices[j], taking into account that `closest' will be inserted into the merged path, thereby remapping the relevant index in the merged path */ new_parent_best_indices[j] = PARENT_VERTEX_IN_MERGED_PATH(nearest_index_in_parent, parent_best_indices[closest], merged_path->num_segments, child_best_indices[closest], parent->children[closest]->num_segments); } } /* do the actual insertion, by adding a pair of lineto's between closest vertices; flag child as inserted */ insert_subpath (merged_path->segments, parent->children[closest]->segments, merged_path->num_segments, parent->children[closest]->num_segments, parent_best_indices[closest], child_best_indices[closest]); merged_path->num_segments += (parent->children[closest]->num_segments + 1); parent->children[closest]->inserted = true; /* update the old arrays to take insertion into account: replace them by the new ones */ for (j = 0; j < parent->num_children; j++) { parent_to_child_distances[j] = new_parent_to_child_distances[j]; parent_best_indices[j] = new_parent_best_indices[j]; child_best_indices[j] = new_child_best_indices[j]; } free (new_parent_to_child_distances); free (new_parent_best_indices); free (new_child_best_indices); } /* End of loop over all children of parent subpath; all >=1 children have now been inserted into the parent, i.e. into the `merged path' which the parent initialized. However, the merged path's segments are all lines; so change the first to a moveto. */ merged_path->segments[0].type = S_MOVETO; merged_paths[i] = merged_path; /* NOTE: SHOULD ALSO REPLACE LAST LINE SEGMENT BY A CLOSEPATH! */ /* delete temporary storage for `closest vertex pairs' and inter-path distances */ free (parent_to_child_distances); free (parent_best_indices); free (child_best_indices); } /* end of loop over parent subpaths */ /* no more annotated paths needed */ delete_subpath_array (annotated_paths, num_paths); return merged_paths; }