/* * Geo.c -- Implementation of common geometric routines. * * Authors : Patrick Lecoanet. * Creation date : * * $Id$ */ /* * Copyright (c) 1993 - 2005 CENA, Patrick Lecoanet -- * * See the file "Copyright" for information on usage and redistribution * of this file, and for a DISCLAIMER OF ALL WARRANTIES. * */ /* * Much of the code here is inspired by (or copied from) the Tk code. */ #include "Geo.h" #include "WidgetInfo.h" #include static const char rcsid[] = "$Id$"; static const char compile_id[]="$Compile: " __FILE__ " " __DATE__ " " __TIME__ " $"; void ZnPolyInit(ZnPoly *poly) { poly->num_contours = 0; poly->contours = NULL; } void ZnPolyContour1(ZnPoly *poly, ZnPoint *pts, unsigned int num_pts, ZnBool cw) { poly->num_contours = 1; poly->contours = &poly->contour1; poly->contour1.num_points = num_pts; poly->contour1.points = pts; poly->contour1.cw = cw; poly->contour1.controls = NULL; } void ZnPolySet(ZnPoly *poly1, ZnPoly *poly2) { ZnPolyFree(poly1); if (poly2->num_contours == 1) { ZnPolyContour1(poly1, poly2->contours[0].points, poly2->contours[0].num_points, poly2->contours[0].cw); if (poly2->contours != &poly2->contour1) { ZnFree(poly2->contours); } } else { poly1->num_contours = poly2->num_contours; poly1->contours = poly2->contours; } } void ZnPolyFree(ZnPoly *poly) { if (poly->num_contours) { unsigned int i; for (i = 0; i < poly->num_contours; i++) { ZnFree(poly->contours[i].points); /* if (poly->contours[i].controls) { ZnFree(poly->contours[i].controls); }*/ } if (poly->contours != &poly->contour1) { ZnFree(poly->contours); } poly->num_contours = 0; poly->contours = NULL; } } void ZnTriStrip1(ZnTriStrip *tristrip, ZnPoint *pts, unsigned int num_pts, ZnBool fan) { tristrip->num_strips = 1; tristrip->strips = &tristrip->strip1; tristrip->strip1.points = pts; tristrip->strip1.num_points = num_pts; tristrip->strip1.fan = fan; } void ZnTriFree(ZnTriStrip *tristrip) { if (tristrip->num_strips) { unsigned int i; for (i = 0; i < tristrip->num_strips; i++) { ZnFree(tristrip->strips[i].points); } if (tristrip->strips != &tristrip->strip1) { ZnFree(tristrip->strips); } tristrip->num_strips = 0; tristrip->strips = NULL; } } /* * Compute the origin of the rectangle given * by position, anchor, width and height. */ void ZnAnchor2Origin(ZnPoint *position, ZnDim width, ZnDim height, Tk_Anchor anchor, ZnPoint *origin) { switch (anchor) { case TK_ANCHOR_CENTER: origin->x = position->x - width/2; origin->y = position->y - height/2; break; case TK_ANCHOR_NW: *origin = *position; break; case TK_ANCHOR_N: origin->x = position->x - width/2; origin->y = position->y; break; case TK_ANCHOR_NE: origin->x = position->x - width; origin->y = position->y; break; case TK_ANCHOR_E: origin->x = position->x - width; origin->y = position->y - height/2; break; case TK_ANCHOR_SE: origin->x = position->x - width; origin->y = position->y - height; break; case TK_ANCHOR_S: origin->x = position->x - width/2; origin->y = position->y - height; break; case TK_ANCHOR_SW: origin->x = position->x; origin->y = position->y - height; break; case TK_ANCHOR_W: origin->x = position->x; origin->y = position->y - height/2; break; } } /* * Compute the anchor position given the bbox origin, width, * height and the anchor. */ void ZnOrigin2Anchor(ZnPoint *origin, ZnDim width, ZnDim height, Tk_Anchor anchor, ZnPoint *position) { switch (anchor) { case TK_ANCHOR_CENTER: position->x = origin->x + width/2; position->y = origin->y + height/2; break; case TK_ANCHOR_NW: *position = *origin; break; case TK_ANCHOR_N: position->x = origin->x + width/2; position->y = origin->y; break; case TK_ANCHOR_NE: position->x = origin->x + width; position->y = origin->y; break; case TK_ANCHOR_E: position->x = origin->x + width; position->y = origin->y + height/2; break; case TK_ANCHOR_SE: position->x = origin->x + width; position->y = origin->y + height; break; case TK_ANCHOR_S: position->x = origin->x + width/2; position->y = origin->y + height; break; case TK_ANCHOR_SW: position->x = origin->x; position->y = origin->y + height; break; case TK_ANCHOR_W: position->x = origin->x; position->y = origin->y + height/2; break; } } /* * Compute the anchor position given a rectangle and * the anchor. The rectangle vertices must be ordered * as for a triangle strip: * * v0 ------------ v2 * | | * | | * v1 ------------ v3 */ void ZnRectOrigin2Anchor(ZnPoint *rect, Tk_Anchor anchor, ZnPoint *position) { switch (anchor) { case TK_ANCHOR_CENTER: position->x = (rect[0].x + rect[3].x) / 2.0; position->y = (rect[0].y + rect[3].y) / 2.0; break; case TK_ANCHOR_NW: *position = *rect; break; case TK_ANCHOR_N: position->x = (rect[0].x + rect[2].x) / 2.0; position->y = (rect[0].y + rect[2].y) / 2.0; break; case TK_ANCHOR_NE: *position = rect[2]; break; case TK_ANCHOR_E: position->x = (rect[2].x + rect[3].x) / 2.0; position->y = (rect[2].y + rect[3].y) / 2.0; break; case TK_ANCHOR_SE: *position = rect[3]; break; case TK_ANCHOR_S: position->x = (rect[1].x + rect[3].x) / 2.0; position->y = (rect[1].y + rect[3].y) / 2.0; break; case TK_ANCHOR_SW: *position = rect[1]; break; case TK_ANCHOR_W: position->x = (rect[0].x + rect[1].x) / 2.0; position->y = (rect[0].y + rect[1].y) / 2.0; break; } } void ZnBBox2XRect(ZnBBox *bbox, XRectangle *r) { r->x = ZnNearestInt(bbox->orig.x); r->y = ZnNearestInt(bbox->orig.y); r->width = ZnNearestInt(bbox->corner.x) - r->x; r->height = ZnNearestInt(bbox->corner.y) - r->y; } void ZnGetStringBBox(char *str, Tk_Font font, ZnPos x, ZnPos y, ZnBBox *str_bbox) { Tk_FontMetrics fm; str_bbox->orig.x = x; str_bbox->corner.x = x + Tk_TextWidth(font, str, (int) strlen(str)); Tk_GetFontMetrics(font, &fm); str_bbox->orig.y = y - fm.ascent; str_bbox->corner.y = str_bbox->orig.y + fm.ascent + fm.descent; } void ZnResetBBox(ZnBBox *bbox) { bbox->orig.x = bbox->orig.y = 0; bbox->corner = bbox->orig; } void ZnCopyBBox(ZnBBox *bbox_from, ZnBBox *bbox_to) { bbox_to->orig = bbox_from->orig; bbox_to->corner = bbox_from->corner; } void ZnIntersectBBox(ZnBBox *bbox1, ZnBBox *bbox2, ZnBBox *bbox_inter) { if ((bbox1->corner.x < bbox2->orig.x) || (bbox1->corner.y < bbox2->orig.y) || (bbox2->corner.x < bbox1->orig.x) || (bbox2->corner.y < bbox1->orig.y)) { ZnResetBBox(bbox_inter); } else { bbox_inter->orig.x = MAX(bbox1->orig.x, bbox2->orig.x); bbox_inter->orig.y = MAX(bbox1->orig.y, bbox2->orig.y); bbox_inter->corner.x = MIN(bbox1->corner.x, bbox2->corner.x); bbox_inter->corner.y = MIN(bbox1->corner.y, bbox2->corner.y); } } ZnBool ZnIsEmptyBBox(ZnBBox *bbox) { return (bbox->orig.x >= bbox->corner.x) || (bbox->orig.y >= bbox->corner.y); } void ZnAddBBoxToBBox(ZnBBox *bbox, ZnBBox *bbox2) { if (ZnIsEmptyBBox(bbox2)) { return; } if (ZnIsEmptyBBox(bbox)) { ZnCopyBBox(bbox2, bbox); } else { bbox->orig.x = MIN(bbox->orig.x, bbox2->orig.x); bbox->orig.y = MIN(bbox->orig.y, bbox2->orig.y); bbox->corner.x = MAX(bbox->corner.x, bbox2->corner.x); bbox->corner.y = MAX(bbox->corner.y, bbox2->corner.y); } } void ZnAddPointToBBox(ZnBBox *bbox, ZnPos px, ZnPos py) { if (ZnIsEmptyBBox(bbox)) { bbox->orig.x = px; bbox->orig.y = py; bbox->corner.x = bbox->orig.x + 1; bbox->corner.y = bbox->orig.y + 1; } else { bbox->orig.x = MIN(bbox->orig.x, px); bbox->orig.y = MIN(bbox->orig.y, py); bbox->corner.x = MAX(bbox->corner.x, px + 1); bbox->corner.y = MAX(bbox->corner.y, py + 1); } } void ZnAddPointsToBBox(ZnBBox *bbox, ZnPoint *points, unsigned int num_points) { ZnReal x1, y1, x2, y2, cur; if (points == NULL) { return; } if (num_points == 0) { return; } if (ZnIsEmptyBBox(bbox)) { x1 = points->x; y1 = points->y; x2 = x1 + 1; y2 = y1 + 1; num_points--; points++; } else { x1 = bbox->orig.x; y1 = bbox->orig.y; x2 = bbox->corner.x; y2 = bbox->corner.y; } for ( ; num_points > 0; num_points--, points++) { cur = points->x; if (cur < x1) { x1 = cur; } if (cur > x2) { x2 = cur; } cur = points->y; if (cur < y1) { y1 = cur; } if (cur > y2) { y2 = cur; } } bbox->orig.x = x1; bbox->orig.y = y1; if (x1 == x2) { x2++; } if (y1 == y2) { y2++; } bbox->corner.x = x2; bbox->corner.y = y2; } void ZnAddStringToBBox(ZnBBox *bbox, char *str, Tk_Font font, ZnPos cx, ZnPos cy) { ZnBBox str_bbox; ZnGetStringBBox(str, font, cx, cy, &str_bbox); ZnAddBBoxToBBox(bbox, &str_bbox); } ZnBool ZnPointInBBox(ZnBBox *bbox, ZnPos x, ZnPos y) { return ((x >= bbox->orig.x) && (x < bbox->corner.x) && (y >= bbox->orig.y) && (y < bbox->corner.y)); } /* * Tell where aa area is with respect to another area. * Return -1 if the first is entirely outside the second, * 1 if it is entirely inside and 0 otherwise. */ int ZnBBoxInBBox(ZnBBox *bbox1, ZnBBox *bbox2) { if ((bbox1->corner.x <= bbox2->orig.x) || (bbox1->orig.x >= bbox2->corner.x) || (bbox1->corner.y <= bbox2->orig.y) || (bbox1->orig.y >= bbox2->corner.y)) { return -1; } if ((bbox2->orig.x <= bbox1->orig.x) && (bbox1->corner.x <= bbox2->corner.x) && (bbox2->orig.y <= bbox1->orig.y) && (bbox1->corner.y <= bbox2->corner.y)) { return 1; } return 0; } /* * Tell where a line is with respect to an area. * Return -1 if the line is entirely outside the bbox, 1 * if it is entirely inside and 0 otherwise. */ int ZnLineInBBox(ZnPoint *p1, ZnPoint *p2, ZnBBox *bbox) { ZnBool p1_inside = ZnPointInBBox(bbox, p1->x, p1->y); ZnBool p2_inside = ZnPointInBBox(bbox, p2->x, p2->y); if (p1_inside != p2_inside) { return 0; } if (p1_inside && p2_inside) { return 1; } /* * Segment may intersect area, check it more thoroughly. */ /* Vertical line */ if (p1->x == p2->x) { if (((p1->y >= bbox->orig.y) ^ (p2->y >= bbox->orig.y)) && (p1->x >= bbox->orig.x) && (p1->x <= bbox->corner.x)) { return 0; } } /* Horizontal line */ else if (p1->y == p2->y) { if (((p1->x >= bbox->orig.x) ^ (p2->x >= bbox->orig.x)) && (p1->y >= bbox->orig.y) && (p1->y <= bbox->corner.y)) { return 0; } } /* Diagonal, do it the hard way. */ else { ZnReal slope = (p2->y - p1->y) / (p2->x - p1->x); ZnDim low, high, x, y; ZnDim bbox_width = bbox->corner.x - bbox->orig.x; ZnDim bbox_height = bbox->corner.y - bbox->orig.y; /* Check against left edge */ if (p1->x < p2->x) { low = p1->x; high = p2->x; } else { low = p2->x; high = p1->x; } y = p1->y + (bbox->orig.x - p1->x) * slope; if ((bbox->orig.x >= low) && (bbox->orig.x <= high) && (y >= bbox->orig.y) && (y <= bbox->corner.y)) return 0; /* Check against right edge */ y += bbox_width * slope; if ((y >= bbox->orig.y) && (y <= bbox->corner.y) && (bbox->corner.x >= low) && (bbox->corner.x <= high)) return 0; /* Check against bottom edge */ if (p1->y < p2->y) { low = p1->y; high = p2->y; } else { low = p2->y; high = p1->y; } x = p1->x + (bbox->orig.y - p1->y) / slope; if ((x >= bbox->orig.x) && (x <= bbox->corner.x) && (bbox->orig.y >= low) && (bbox->orig.y <= high)) return 0; /* Check against top edge */ x += bbox_height / slope; if ((x >= bbox->orig.x) && (x <= bbox->corner.x) && (bbox->corner.y >= low) && (bbox->corner.y <= high)) return 0; } return -1; } ZnBool ZnTestCCW(ZnPoint *points, unsigned int num_points) { ZnPoint *p, *p_p=NULL, *p_n=NULL, min; ZnReal xprod; unsigned int i, min_index; if (num_points < 3) { return True; } /* * First find the lowest rightmost vertex. In X11 this is the * topmost one. */ p = points; min = *p; min_index = 0; for (i = 1, p++; i < num_points; i++, p++) { if ((p->y < min.y) || ((p->y == min.y) && (p->x > min.x))) { min_index = i; min = *p; } } /* * Then find the indices of the previous and next * vertices. */ p = &points[min_index]; /*printf("min index %d, prev %d, next %d\n", min_index, (min_index+(num_points-1))%num_points, (min_index+1)%num_points);*/ /*printf("lower point index %d %d %d\n", (min_index+(num_points-1))%num_points, min_index, (min_index+1)%num_points);*/ /* * Try to find preceding and following points that are not * the same as the base point. */ for (i = 1; i < num_points; i++) { p_p = &points[(min_index+(num_points-i))%num_points]; /* min_index-1 */ if ((p_p->x != p->x) || (p_p->y != p->y)) { break; } } for (i = 1; i < num_points; i++) { p_n = &points[(min_index+i)%num_points]; if ((p_p->x != p->x) || (p_p->y != p->y)) { break; } } xprod = ((p->x - p_p->x) * (p_n->y - p->y)) - ((p->y - p_p->y) * (p_n->x - p->x)); return (xprod < 0.0); /* Should be > 0 but X11 has Y axis reverted. */ } /* * ZnShiftLine -- * Given two points describing a line and a distance, return * to points describing a line parallel to it at the given distance. * When looking the line from p1 to p2 the new line will be dist away * on its left. Negative values are allowed for dist, resulting in a line * on the right. */ void ZnShiftLine(ZnPoint *p1, ZnPoint *p2, ZnReal dist, ZnPoint *p3, ZnPoint *p4) { static int shift_table[129]; ZnBool dx_neg, dy_neg; int dx, dy; /* * Initialize the conversion table. */ if (shift_table[0] == 0) { int i; ZnReal tangent, cosine; for (i = 0; i <= 128; i++) { tangent = i/128.0; cosine = 128/cos(atan(tangent)) + 0.5; shift_table[i] = (int) cosine; } } *p3 = *p1; dx = (int) (p2->x - p1->x); dy = (int) (p2->y - p1->y); if (dx < 0) { dx = -dx; dx_neg = True; } else { dx_neg = False; } if (dy < 0) { dy = -dy; dy_neg = True; } else { dy_neg = False; } if ((dy < PRECISION_LIMIT) && (dx < PRECISION_LIMIT)) { fprintf(stderr, "ShiftLine: segment is a point\n"); return; } if (dy <= dx) { dy = (((int) dist * shift_table[(dy*128)/dx]) + 64) / 128; if (!dx_neg) { dy = -dy; } p3->y += dy; } else { dx = (((int) dist * shift_table[(dx*128)/dy]) + 64) / 128; if (dy_neg) { dx = -dx; } p3->x += dx; } p4->x = p3->x + (p2->x - p1->x); p4->y = p3->y + (p2->y - p1->y); } /* * IntersectLines -- * Given two lines described by two points, compute their intersection. * The function returns True if the lines are not parallel and False * otherwise. */ ZnBool ZnIntersectLines(ZnPoint *a1, ZnPoint *a2, ZnPoint *b1, ZnPoint *b2, ZnPoint *pi) { ZnReal dxadyb, dxbdya, dxadxb, dyadyb, p, q; dxadyb = (a2->x - a1->x)*(b2->y - b1->y); dxbdya = (b2->x - b1->x)*(a2->y - a1->y); dxadxb = (a2->x - a1->x)*(b2->x - b1->x); dyadyb = (a2->y - a1->y)*(b2->y - b1->y); if (dxadyb == dxbdya) { return False; } p = a1->x*dxbdya - b1->x*dxadyb + (b1->y - a1->y)*dxadxb; q = dxbdya - dxadyb; if (q < 0) { p = -p; q = -q; } if (p < 0) { pi->x = - ((-p + q/2)/q); } else { pi->x = (p + q/2)/q; } p = a1->y*dxadyb - b1->y*dxbdya + (b1->x - a1->x)*dyadyb; q = dxadyb - dxbdya; if (q < 0) { p = -p; q = -q; } if (p < 0) { pi->y = - ((-p + q/2)/q); } else { pi->y = (p + q/2)/q; } return True; } /* * InsetPolygon -- * Inset the given polygon by the given amount. The * value can be negative, in this case the polygon will * be outset. */ /**** A FINIR ****/ void ZnInsetPolygon(ZnPoint *p, unsigned int num_points, ZnDim inset) { ZnPoint *p1, *p2; ZnPoint new_p1, new_p2; /* ZnPoint shift1, shift2;*/ unsigned int i, processed_points; processed_points = 0; if ((p->x == p[num_points-1].x) && (p->y == p[num_points-1].y)) { num_points--; } for (p1 = p, p2 = p1+1, i = 0; i < num_points; i++, p1 = p2, p2++) { /* * Wrap to the first point. */ if (i == num_points-1) { p2 = p; } /* * Skip duplicate vertices. */ if ((p2->x == p1->x) && (p2->y == p1->y)) { continue; } ZnShiftLine(p1, p2, inset, &new_p1, &new_p2); if (processed_points >= 1) { } } } /* * Compute the two corner points of a thick line end. * Two points describing the line segment and the width * are given as input. If projecting is true this function * mimics the X11 line projecting behaviour. The computed * end is located around p2. */ void ZnGetButtPoints(ZnPoint *p1, ZnPoint *p2, ZnDim width, ZnBool projecting, ZnPoint *c1, ZnPoint *c2) { ZnReal w_2 = width/2.0; ZnDim length = hypot(p2->x - p1->x, p2->y - p1->y); ZnReal delta_x, delta_y; if (length == 0.0) { c1->x = c2->x = p2->x; c1->y = c2->y = p2->y; } else { delta_x = -w_2 * (p2->y - p1->y) / length; delta_y = w_2 * (p2->x - p1->x) / length; c1->x = p2->x + delta_x; c2->x = p2->x - delta_x; c1->y = p2->y + delta_y; c2->y = p2->y - delta_y; if (projecting) { c1->x += delta_y; c2->x += delta_y; c1->y -= delta_x; c2->y -= delta_x; } } } /* * Compute the inside and outside points of the mitered * corner formed by a thick line going through 3 points. * If the angle formed by the three points is less than * 11 degrees, False is returned an no points are computed. * Else True is returned and the points are in c1, c2. * * If someday the code is switched to REAL coordinates, we * must round each coordinate to the nearer integer to mimic * the way pixels are drawn. Sample code: floor(p->x+0.5); * * Hmmm, the switch has been done but not the rounding ;-) */ ZnBool ZnGetMiterPoints(ZnPoint *p1, ZnPoint *p2, ZnPoint *p3, ZnDim width, ZnPoint *c1, ZnPoint *c2) { static ZnReal deg11 = (11.0*2.0*M_PI)/360.0; ZnReal theta1; /* angle of p2-p1 segment. */ ZnReal theta2; /* angle of p2-p3 segment. */ ZnReal theta; /* angle of the joint */ ZnReal theta3; /* angle of bisector of the joint toward * the external point of the joint. */ ZnReal dist; /* distance of the external points * of the corner from the mid point * p2. */ ZnReal delta_x, delta_y; /* projection of (dist,theta3) on x * and y. */ if (p2->y == p1->y) { theta1 = (p2->x < p1->x) ? 0.0 : M_PI; } else if (p2->x == p1->x) { theta1 = (p2->y < p1->y) ? M_PI/2.0 : -M_PI/2.0; } else { theta1 = atan2(p1->y - p2->y, p1->x - p2->x); } if (p3->y == p2->y) { theta2 = (p3->x > p2->x) ? 0.0 : M_PI; } else if (p3->x == p2->x) { theta2 = (p3->y > p2->y) ? M_PI/2.0 : -M_PI/2.0; } else { theta2 = atan2(p3->y - p2->y, p3->x - p2->x); } theta = theta1 - theta2; if (theta > M_PI) { theta -= 2.0*M_PI; } else if (theta < -M_PI) { theta += 2*M_PI; } if ((theta < deg11) && (theta > -deg11)) { return False; } /* * Compute the distance of the internal and external * corner points from the intersection p2 (considered * at 0,0). */ dist = 0.5*width / sin(0.5*theta); dist = ABS(dist); /* * Compute the angle of the bisector of the joint that * goes toward the outside of the joint (the left hand * when looking from p1-p2). */ theta3 = (theta1 + theta2)/2.0; if (sin(theta3 - (theta1 + M_PI)) < 0.0) { theta3 += M_PI; } delta_x = dist * cos(theta3); c1->x = p2->x + delta_x; c2->x = p2->x - delta_x; delta_y = dist * sin(theta3); c1->y = p2->y + delta_y; c2->y = p2->y - delta_y; return True; } /* * Tell where a thick polyline is with respect to an area. * Return -1 if the polyline is entirely outside the bbox, 1 * if it is entirely inside and 0 otherwise. The joints can * be specified as JoinMiter, JoinRound, JoinBevel. The ends * can be: CapRound, CapButt, CapProjecting. */ int ZnPolylineInBBox(ZnPoint *points, unsigned int num_points, ZnDim width, int cap_style, int join_style, ZnBBox *bbox) { unsigned int count; int inside = -1; ZnBool do_miter_as_bevel; ZnPoint poly[4]; /* * If the first point is inside the area, change inside * accordingly. */ if ((points[0].x >= bbox->orig.x) && (points[0].x <= bbox->corner.x) && (points[0].y >= bbox->orig.y) && (points[0].y <= bbox->corner.y)) { inside = 1; } /* * Now iterate through all the edges. Compute a polygon for * each and test it against the area. At each vertex an oval * of radius width/2 is also tested to account for round ends * and joints. */ do_miter_as_bevel = False; for (count = num_points; count >= 2; count--, points++) { /* * Test a circle around the first point if CapRound or * around every joint if JoinRound. */ if (((cap_style == CapRound) && (count == num_points)) || ((join_style == JoinRound) && (count != num_points))) { if (ZnOvalInBBox(points, width, width, bbox) != inside) { return 0; } } /* * Build a polygon to represent an edge from the current * point to the next. Special cases for the first and * last edges to implement the line ends. */ /* * First vertex of the edge */ if (count == num_points) { ZnGetButtPoints(&points[1], points, width, cap_style == CapProjecting, poly, &poly[1]); } /* * Here we are at a joint starting a new edge. If the * joint is mitered, start by carrying over the points * from the previous edge. Otherwise compute new points * for a butt end. */ else if ((join_style == JoinMiter) && !do_miter_as_bevel) { poly[0] = poly[3]; poly[1] = poly[2]; } else { ZnGetButtPoints(&points[1], points, width, 0, poly, &poly[1]); /* * If the previous joint was beveled (or considered so), * check the polygon that fill the bevel. It has more or * less an X shape, i.e, it's self intersecting. If this * is not ok, it may be necessary to permutte poly[1] & * poly[2]. */ if ((join_style == JoinBevel) || do_miter_as_bevel) { if (ZnPolygonInBBox(poly, 4, bbox, NULL) != inside) { return 0; } do_miter_as_bevel = False; } } /* * Opposite vertex of the edge. */ if (count == 2) { ZnGetButtPoints(points, &points[1], width, cap_style == CapProjecting, &poly[2], &poly[3]); } else if (join_style == JoinMiter) { if (ZnGetMiterPoints(points, &points[1], &points[2], width, &poly[2], &poly[3]) == False) { do_miter_as_bevel = True; ZnGetButtPoints(points, &points[1], width, 0, &poly[2], &poly[3]); } } else { ZnGetButtPoints(points, &points[1], width, 0, &poly[2], &poly[3]); } if (ZnPolygonInBBox(poly, 4, bbox, NULL) != inside) { return 0; } } /* * Test a circle around the last point if CapRound. */ if (cap_style == CapRound) { if (ZnOvalInBBox(points, width, width, bbox) != inside) { return 0; } } return inside; } /* * Tell where a polygon is with respect to an area. * Return -1 if the polygon is entirely outside the bbox, 1 * if it is entirely inside and 0 otherwise. If area_enclosed * is not NULL it tells whether the area is enclosed by the * polygon or not. */ int ZnPolygonInBBox(ZnPoint *points, unsigned int num_points, ZnBBox *bbox, ZnBool *area_enclosed) { int inside, count; ZnPoint *p, *head, *first, *second; ZnBool closed; if (area_enclosed) { *area_enclosed = False; } p = head = points; closed = True; count = num_points-2; /* * Check to see if closed. If not adjust num_points and * record this. */ if ((points[0].x != points[num_points-1].x) || (points[0].y != points[num_points-1].y)) { closed = False; count = num_points-1; } /* * Get the status of the first edge. */ inside = ZnLineInBBox(&p[0], &p[1], bbox); p++; if (inside == 0) { return 0; } for (; count > 0; p++, count--) { first = &p[0]; /* * Pretend the polygon is closed if this is not the case. */ if (!closed && (count == 1)) { second = head; } else { second = &p[1]; } if (ZnLineInBBox(first, second, bbox) != inside) { return 0; } } /* * If all the edges are inside the area, the polygon is * inside the area. If all the edges are outside, the polygon * may completely enclose the area. Test if the origin of * the area is inside the polygon to decide this. */ if (inside == 1) { return 1; } /*printf("PolygonInBBox, np = %d, x = %g, y = %g, dist = %g\n", num_points, bbox->orig.x, bbox->orig.y, PolygonToPointDist(points, num_points, &bbox->orig));*/ if (ZnPolygonToPointDist(points, num_points, &bbox->orig) <= 0.0) { if (area_enclosed) { *area_enclosed = True; } return 0; } return -1; } /* * Tell where an oval is with respect to an area. * Return -1 if the oval is entirely outside the bbox, 1 * if it is entirely inside and 0 otherwise. */ int ZnOvalInBBox(ZnPoint *center, ZnDim width, ZnDim height, ZnBBox *bbox) { ZnPoint origin, corner; ZnDim w_2, h_2; ZnReal delta_x, delta_y; w_2 = (width+1)/2; h_2 = (height+1)/2; origin.x = center->x - w_2; origin.y = center->y - h_2; corner.x = center->x + w_2; corner.y = center->y + h_2; /* * This if the oval bbox is completely inside or outside * of the area. Trivial case first. */ if ((bbox->orig.x <= origin.x) && (bbox->corner.x >= corner.x) && (bbox->orig.y <= origin.y) && (bbox->corner.y >= corner.y)) { return 1; } if ((bbox->corner.x < origin.x) || (bbox->orig.x > corner.x) || (bbox->corner.y < origin.y) || (bbox->orig.y > corner.y)) { return -1; } /* * Then test all sides of the area against the oval center. * If the point of a side closest to the center is inside * the oval, then the oval intersects the area. */ /* * Compute the square of the Y axis distance, reducing * the oval to a unit circle to take into account the * shape factor. */ delta_y = bbox->orig.y - center->y; if (delta_y < 0.0) { delta_y = center->y - bbox->corner.y; if (delta_y < 0.0) { delta_y = 0.0; } } delta_y /= h_2; delta_y *= delta_y; /* * Test left and then right edges. */ delta_x = (bbox->orig.x - center->x) / w_2; delta_x *= delta_x; if ((delta_x + delta_y) <= 1.0) { return 0; } delta_x = (bbox->corner.x - center->x) / w_2; delta_x *= delta_x; if ((delta_x + delta_y) <= 1.0) { return 0; } /* * Compute the square of the X axis distance, reducing * the oval to a unit circle to take into account the * shape factor. */ delta_x = bbox->orig.x - center->x; if (delta_x < 0.0) { delta_x = center->x - bbox->corner.x; if (delta_x < 0.0) { delta_x = 0.0; } } delta_x /= w_2; delta_x *= delta_x; /* * Test top and then bottom edges. */ delta_y = (bbox->orig.y - center->y) / h_2; delta_y *= delta_y; if ((delta_x + delta_y) <= 1.0) { return 0; } delta_y = (bbox->corner.y - center->y) / h_2; delta_y *= delta_y; if ((delta_x + delta_y) <= 1.0) { return 0; } return -1; } /* * Tell if a point is in an angular range whose center is 0,0. * The range is specified by a starting angle and an angle extent. * The use of a double here is important, don't change it. In some * case we need to normalize a figure to take care of its shape and * the result needs precision. */ ZnBool ZnPointInAngle(int start_angle, int angle_extent, ZnPoint *p) { ZnReal point_angle; int angle_diff; if ((p->x == 0) && (p->y == 0)) { point_angle = 0.0; } else { point_angle = atan2(p->y, p->x) * 180.0 / M_PI; } angle_diff = (ZnNearestInt(point_angle) - start_angle) % 360; if (angle_diff < 0) { angle_diff += 360; } return ((angle_diff <= angle_extent) || ((angle_extent < 0) && ((angle_diff - 360) >= angle_extent))); } /* * PointPolarToCartesian -- * Convert a point in polar coordinates (rho, theta) * in a reference system described by angle heading * (angles running clockwise) to a point in cartesian * coordinates (delta_x, delta_y). * */ void ZnPointPolarToCartesian(ZnReal heading, ZnReal rho, ZnReal theta, ZnReal *delta_x, ZnReal *delta_y) { ZnReal to_angle; /* Compute angle in trigonometric system */ /* to_angle = ZnDegRad(theta) + heading - M_PI_2;*/ to_angle = heading - ZnDegRad(theta) - M_PI_2; /* Compute cartesian coordinates */ *delta_x = rho * cos(to_angle); *delta_y = rho * sin(to_angle); } /* * Return a vector angle given its projections */ ZnReal ZnProjectionToAngle(ZnReal dx, ZnReal dy) { if (dx == 0) { if (dy < 0) { return -M_PI_2; } else if (dy > 0) { return M_PI_2; } else { return 0.0; } } else if (dx < 0) { return atan(dy / dx) - M_PI; } else { return atan(dy / dx); } return 0.0; } /* * Tell if an horizontal line intersect an axis aligned * elliptical arc. * * Returns True if the line described by (x1, x2, y) intersects * the arc described by (r1, r2, start_angle and angle_extent). * This arc is origin centered. */ ZnBool ZnHorizLineToArc(ZnReal x1, ZnReal x2, ZnReal y, ZnReal rx, ZnReal ry, int start_angle, int angle_extent) { ZnReal tmp, x; ZnPoint t; /* * Compute the x-coordinate of one possible intersection point * between the arc and the line. Use a transformed coordinate * system where the oval is a unit circle centered at the origin. * Then scale back to get actual x-coordinate. */ t.y = y/ry; tmp = 1 - t.y*t.y; if (tmp < 0.0) { return False; } t.x = sqrt(tmp); x = t.x*rx; /* * Test both intersection points. */ if ((x >= x1) && (x <= x2) && ZnPointInAngle((int) start_angle, (int) angle_extent, &t)) { return True; } t.x = -t.x; if ((-x >= x1) && (-x <= x2) && ZnPointInAngle((int) start_angle, (int) angle_extent, &t)) { return True; } return False; } /* * Tell if an vertical line intersect an axis aligned * elliptical arc. * * Returns True if the line described by (x1, x2, y) intersects * the arc described by (r1, r2, start_angle and angle_extent). * This arc is origin centered. */ ZnBool ZnVertLineToArc(ZnReal x, ZnReal y1, ZnReal y2, ZnReal rx, ZnReal ry, int start_angle, int angle_extent) { ZnReal tmp, y; ZnPoint t; /* * Compute the y-coordinate of one possible intersection point * between the arc and the line. Use a transformed coordinate * system where the oval is a unit circle centered at the origin. * Then scale back to get actual y-coordinate. */ t.x = x/rx; tmp = 1 - t.x*t.x; if (tmp < 0.0) { return False; } t.y = sqrt(tmp); y = t.y*ry; /* * Test both intersection points. */ if ((y > y1) && (y < y2) && ZnPointInAngle((int) start_angle, (int) angle_extent, &t)) { return True; } t.y = -t.y; if ((-y > y1) && (-y < y2) && ZnPointInAngle((int) start_angle, (int) angle_extent, &t)) { return True; } return False; } /* * Return the distance of the given point to the rectangle * described by rect. Return negative values for points in * the rectangle. */ ZnDim ZnRectangleToPointDist(ZnBBox *bbox, ZnPoint *p) { ZnDim new_dist, dist; ZnPoint p1, p2; p1.x = bbox->orig.x; p1.y = p2.y = bbox->orig.y; p2.x = bbox->corner.x; dist = ZnLineToPointDist(&p1, &p2, p, NULL); if (dist == 0.0) { return 0.0; } p1 = p2; p2.y = bbox->corner.y; new_dist = ZnLineToPointDist(&p1, &p2, p, NULL); dist = MIN(dist, new_dist); if (dist == 0.0) { return 0.0; } p1 = p2; p2.x = bbox->orig.x; new_dist = ZnLineToPointDist(&p1, &p2, p, NULL); dist = MIN(dist, new_dist); if (dist == 0.0) { return 0.0; } p1 = p2; p2.y = bbox->orig.y; new_dist = ZnLineToPointDist(&p1, &p2, p, NULL); dist = MIN(dist, new_dist); if (ZnPointInBBox(bbox, p->x, p->y)) { return -dist; } else { return dist; } } /* * Return the distance of the given point to the line * described by , .. */ ZnDim ZnLineToPointDist(ZnPoint *p1, ZnPoint *p2, ZnPoint *p, ZnPoint *closest) { ZnReal x, y; ZnReal x_int, y_int; /* * First compute the closest point on the line. This is done * separately for vertical, horizontal, other lines. */ /* Vertical */ if (p1->x == p2->x) { x = p1->x; if (p1->y >= p2->y) { y_int = MIN(p1->y, p->y); y_int = MAX(y_int, p2->y); } else { y_int = MIN(p2->y, p->y); y_int = MAX(y_int, p1->y); } y = y_int; } /* Horizontal */ else if (p1->y == p2->y) { y = p1->y; if (p1->x >= p2->x) { x_int = MIN(p1->x, p->x); x_int = MAX(x_int, p2->x); } else { x_int = MIN(p2->x, p->x); x_int = MAX(x_int, p1->x); } x = x_int; } /* * Other. Compute its parameters of form y = a1*x + b1 and * then compute the parameters of the perpendicular passing * through the point y = a2*x + b2, last find the closest point * on the segment. */ else { ZnReal a1, a2, b1, b2; a1 = (p2->y - p1->y) / (p2->x - p1->x); b1 = p1->y - a1*p1->x; a2 = -1.0/a1; b2 = p->y - a2*p->x; x = (b2 - b1) / (a1 - a2); y = a1*x + b1; if (p1->x > p2->x) { if (x > p1->x) { x = p1->x; y = p1->y; } else if (x < p2->x) { x = p2->x; y = p2->y; } } else { if (x > p2->x) { x = p2->x; y = p2->y; } else if (x < p1->x) { x = p1->x; y = p1->y; } } } if (closest) { closest->x = x; closest->y = y; } /* Return the distance */ return hypot(p->x - x, p->y - y); } /* * Return the distance of the polygon described by * points, to the given point. If the point is * inside return values are negative. */ ZnDim ZnPolygonToPointDist(ZnPoint *points, unsigned int num_points, ZnPoint *p) { ZnDim best_distance, dist; int intersections; int x_int, y_int; ZnPoint *first_point; ZnReal x, y; ZnPoint p1, p2; /* * The algorithm iterates through all the edges of the polygon * computing for each the distance to the point and whether a vertical * ray starting at the point, intersects the edge. The smallest * distance of all edges is stored in best_distance while intersections * hold the count of edges to rays intersections. For more informations * on how the distance is computed see LineToPointDist. */ best_distance = 1.0e40; intersections = 0; first_point = points; /* * Check to see if closed. Adjust num_points to open it (the * algorithm always consider a set of points as a closed polygon). */ if ((points[0].x == points[num_points-1].x) && (points[0].y == points[num_points-1].y)) { num_points--; } for ( ; num_points >= 1; num_points--, points++) { p1 = points[0]; /* * Wrap over to the first point. */ if (num_points == 1) { p2 = *first_point; } else { p2 = points[1]; } /* * First try to find the closest point on this edge. */ /* Vertical edge */ if (p1.x == p2.x) { x = p1.x; if (p1.y >= p2.y) { y_int = (int) MIN(p1.y, p->y); y_int = (int) MAX(y_int, p2.y); } else { y_int = (int) MIN(p2.y, p->y); y_int = (int) MAX(y_int, p1.y); } y = y_int; } /* Horizontal edge */ else if (p1.y == p2.y) { y = p1.y; if (p1.x >= p2.x) { x_int = (int) MIN(p1.x, p->x); x_int = (int) MAX(x_int, p2.x); if ((p->y < y) && (p->x < p1.x) && (p->x >= p2.x)) { intersections++; } } else { x_int = (int) MIN(p2.x, p->x); x_int = (int) MAX(x_int, p1.x); if ((p->y < y) && (p->x < p2.x) && (p->x >= p1.x)) { intersections++; } } x = x_int; } /* Other */ else { ZnReal a1, b1, a2, b2; a1 = (p2.y - p1.y) / (p2.x - p1.x); b1 = p1.y - a1 * p1.x; a2 = -1.0/a1; b2 = p->y - a2 * p->x; x = (b2 - b1)/(a1 - a2); y = a1 * x + b1; if (p1.x > p2.x) { if (x > p1.x) { x = p1.x; y = p1.y; } else if (x < p2.x) { x = p2.x; y = p2.y; } } else { if (x > p2.x) { x = p2.x; y = p2.y; } else if (x < p1.x) { x = p1.x; y = p1.y; } } if (((a1 * p->x + b1) > p->y) && /* True if point is lower */ (p->x >= MIN(p1.x, p2.x)) && (p->x < MAX(p1.x, p2.x))) { intersections++; } } /* * Now compute the distance to the closest point and * keep it if it is the shortest. */ dist = hypot(p->x - x, p->y - y); best_distance = MIN(best_distance, dist); /* * We can safely escape here if distance is zero. */ if (best_distance == 0.0) { return 0.0; } } /* * Well, all the edges are processed, if the * intersection count is odd the point is inside. */ if (intersections & 0x1) { return -best_distance; } else { return best_distance; } } /* * Return the distance of a thick polyline to the * given point. Cap and Join parameters are considered * in the process. */ ZnDim ZnPolylineToPointDist(ZnPoint *points, unsigned int num_points, ZnDim width, int cap_style, int join_style, ZnPoint *p) { ZnBool miter2bevel = False; unsigned int count; ZnPoint *ptr; ZnPoint outline[5]; ZnDim dist, best_dist, h_width; best_dist = 1.0e36; h_width = width/2.0; for (count = num_points, ptr = points; count >= 2; count--, ptr++) { if (((cap_style == CapRound) && (count == num_points)) || ((join_style == JoinRound) && (count != num_points))) { dist = hypot(ptr->x - p->x, ptr->y - p->y) - h_width; if (dist <= 0.0) { best_dist = 0.0; goto done; } else if (dist < best_dist) { best_dist = dist; } } /* * Build the polygonal outline of the current edge. */ if (count == num_points) { ZnGetButtPoints(&ptr[1], ptr, width, cap_style==CapProjecting, outline, &outline[1]); } else if ((join_style == JoinMiter) && !miter2bevel) { outline[0] = outline[3]; outline[1] = outline[2]; } else { ZnGetButtPoints(&ptr[1], ptr, width, 0, outline, &outline[1]); /* * If joints are beveled, check the distance to the polygon * that fills the joint. */ if ((join_style == JoinBevel) || miter2bevel) { outline[4] = outline[0]; dist = ZnPolygonToPointDist(outline, 5, p); if (dist <= 0.0) { best_dist = 0.0; goto done; } else if (dist < best_dist) { best_dist = dist; } miter2bevel = False; } } if (count == 2) { ZnGetButtPoints(ptr, &ptr[1], width, cap_style==CapProjecting, &outline[2], &outline[3]); } else if (join_style == JoinMiter) { if (ZnGetMiterPoints(ptr, &ptr[1], &ptr[2], width, &outline[2], &outline[3]) == False) { miter2bevel = True; ZnGetButtPoints(ptr, &ptr[1], width, 0, &outline[2], &outline[3]); } /*printf("2=%g+%g, 3=%g+%g\n", outline[2].x, outline[2].y, outline[3].x, outline[3].y);*/ } else { ZnGetButtPoints(ptr, &ptr[1], width, 0, &outline[2], &outline[3]); } outline[4] = outline[0]; /*printf("0=%g+%g, 1=%g+%g, 2=%g+%g, 3=%g+%g, 4=%g+%g\n", outline[0].x, outline[0].y, outline[1].x, outline[1].y, outline[2].x, outline[2].y, outline[3].x, outline[3].y, outline[4].x, outline[4].y);*/ dist = ZnPolygonToPointDist(outline, 5, p); if (dist <= 0.0) { best_dist = 0.0; goto done; } else if (dist < best_dist) { best_dist = dist; } } /* * Test the final point if cap style is round. The code so far * has only handled the butt and projecting cases. */ if (cap_style == CapRound) { dist = hypot(ptr->x - p->x, ptr->y - p->y) - h_width; if (dist <= 0.0) { best_dist = 0.0; goto done; } else if (dist < best_dist) { best_dist = dist; } } done: return best_dist; } /* * Return the distance of the given oval to the point given. * The oval is described by its bounding box , * the thickness of its outline . Return values are negative * if the point is inside. */ ZnDim ZnOvalToPointDist(ZnPoint *center, ZnDim width, ZnDim height, ZnDim line_width, ZnPoint *p) { ZnReal x_delta, y_delta; /* ZnReal x_diameter, y_diameter;*/ ZnDim scaled_distance; ZnDim distance_to_outline; ZnDim distance_to_center; /* * Compute the distance from the point given to the center * of the oval. Then compute the same distance in a coordinate * system where the oval is a circle with unit radius. */ x_delta = p->x - center->x; y_delta = p->y - center->y; distance_to_center = hypot(x_delta, y_delta); scaled_distance = hypot(x_delta / ((width + line_width) / 2.0), y_delta / ((height + line_width) / 2.0)); /* * If the scaled distance is greater than 1.0 the point is outside * the oval. Compute the distance to the edge and convert it back * to the original coordinate system. This distance is not much * accurate and can overestimate the real distance if the oval is * very eccentric. */ if (scaled_distance > 1.0) { distance_to_outline = (distance_to_center / scaled_distance) * (scaled_distance - 1.0); return distance_to_outline; } /* * The point is inside the oval. Compute the distance as above and check * if the point is within the outline. */ if (scaled_distance > 1.0e-10) { distance_to_outline = (distance_to_center / scaled_distance) * (1.0 - scaled_distance) - line_width; } else { /* * If the point is very close to the center avoid dividing by a * very small number, take another method. */ if (width < height) distance_to_outline = (width - line_width) / 2; else distance_to_outline = (height - line_width) / 2; } if (distance_to_outline < 0.0) return 0.0; else return -distance_to_outline; } static int bezier_basis[][4] = { { -1, 3, -3, 1}, { 3, -6, 3, 0}, { -3, 3, 0, 0}, { 1, 0, 0, 0} }; /* ********************************************************************************** * * Arc2Param -- * * Given a Bezier curve describing an arc and an angle return the parameter * value for the intersection point between the arc and the ray at angle. * ********************************************************************************** */ #define EVAL(coeff, t) (((coeff[0]*t + coeff[1])*t + coeff[2]) * t + coeff[3]) static ZnReal Arc2Param(ZnPoint *controls, ZnReal angle) { ZnReal coeff_x[4], coeff_y[4]; ZnReal min_angle, min_t, max_angle, max_t, cur_angle, cur_t; int i, j, depth = 0; /* assume angle >= 0 */ while (angle > M_PI) { angle -= 2 * M_PI; } for (i = 0; i < 4; i++) { coeff_x[i] = coeff_y[i] = 0.0; for (j = 0; j < 4; j++) { coeff_x[i] += bezier_basis[i][j] * controls[j].x; coeff_y[i] += bezier_basis[i][j] * controls[j].y; } } min_angle = atan2(controls[0].y, controls[0].x); max_angle = atan2(controls[3].y, controls[3].x); if (max_angle < min_angle) { min_angle -= M_PI+M_PI; } if (angle > max_angle) { angle -= M_PI+M_PI; } min_t = 0.0; max_t = 1.0; while (depth < 15) { cur_t = (max_t + min_t) / 2.0; cur_angle = atan2(EVAL(coeff_y, cur_t), EVAL(coeff_x, cur_t)); if (angle > cur_angle) { min_t = cur_t; min_angle = cur_angle; } else { max_t = cur_t; max_angle = cur_angle; } depth += 1; } if ((max_angle-angle) < (angle-min_angle)) { return max_t; } return min_t; } #undef EVAL /* ********************************************************************************** * * BezierSubdivide -- * * Subdivide a Bezier curve given by controls at parameter t. Return * in controls, the first or the last part depending on boolean first. * ********************************************************************************** */ static void BezierSubdivide(ZnPoint *controls, ZnReal t, ZnBool first) { ZnReal s = 1.0 - t; ZnPoint r[7]; ZnPoint a; r[0] = controls[0]; r[6] = controls[3]; a.x = s*controls[1].x + t*controls[2].x; a.y = s*controls[1].y + t*controls[2].y; r[1].x = s*r[0].x + t*controls[1].x; r[1].y = s*r[0].y + t*controls[1].y; r[2].x = s*r[1].x + t*a.x; r[2].y = s*r[1].y + t*a.y; r[5].x = s*controls[2].x + t*r[6].x; r[5].y = s*controls[2].y + t*r[6].y; r[4].x = s*a.x + t*r[5].x; r[4].y = s*a.y + t*r[5].y; r[3].x = s*r[2].x + t*r[4].x; r[3].y = s*r[2].y + t*r[4].y; if (first) { memcpy(controls, r, 4*sizeof(ZnPoint)); } else { memcpy(controls, &r[3], 4*sizeof(ZnPoint)); } } /* ********************************************************************************** * * ZnGetBezierPoints -- * Use recursive subdivision to approximate the curve. The subdivision stops * when the error is under eps. * This algorithm is adaptive, meaning that it computes the minimum number * of segments needed to render each curve part. * ********************************************************************************** */ void ZnGetBezierPoints(ZnPoint *p1, ZnPoint *c1, ZnPoint *c2, ZnPoint *p2, ZnList to_points, ZnReal eps) { ZnReal dist; dist = ZnLineToPointDist(p1, p2, c1, NULL); if ((dist < eps) && ((c1->x != c2->x) || (c1->y != c2->y))) { dist = ZnLineToPointDist(p1, p2, c2, NULL); } if (dist > eps) { ZnPoint mid_segm, new_c1, new_c2; /* * Subdivide the curve at t = 0.5 * and compute each new curve. */ mid_segm.x = (p1->x + 3*c1->x + 3*c2->x + p2->x) / 8.0; mid_segm.y = (p1->y + 3*c1->y + 3*c2->y + p2->y) / 8.0; new_c1.x = (p1->x + c1->x) / 2.0; new_c1.y = (p1->y + c1->y) / 2.0; new_c2.x = (p1->x + 2*c1->x + c2->x) / 4.0; new_c2.y = (p1->y + 2*c1->y + c2->y) / 4.0; ZnGetBezierPoints(p1, &new_c1, &new_c2, &mid_segm, to_points, eps); new_c1.x = (c1->x + 2*c2->x + p2->x) / 4.0; new_c1.y = (c1->y + 2*c2->y + p2->y) / 4.0; new_c2.x = (c2->x + (p2->x)) / 2.0; new_c2.y = (c2->y + (p2->y)) / 2.0; ZnGetBezierPoints(&mid_segm, &new_c1, &new_c2, p2, to_points, eps); } else { /* * Flat enough add the end to the current path. * The start should already be there. */ ZnListAdd(to_points, p2, ZnListTail); } } /* ********************************************************************************** * * ZnGetBezierPath -- * Compute in to_points a new set of points describing a Bezier path based * on the control points given in from_points. * If more than four points are given, the algorithm iterate over the * set using the last point of a segment as the first point of the next. * If 3 points are left, they are interpreted as a Bezier segment with * coincident internal control points. If 2 points are left a straight * is emitted. * ********************************************************************************** */ void ZnGetBezierPath(ZnList from_points, ZnList to_points) { ZnPoint *fp; int num_fp, i; fp = ZnListArray(from_points); num_fp = ZnListSize(from_points); /* * make sure the output vector is empty, then add the first point. */ ZnListEmpty(to_points); ZnListAdd(to_points, fp, ZnListTail); for (i = 0; i < num_fp; ) { if (i < (num_fp-3)) { ZnGetBezierPoints(fp, fp+1, fp+2, fp+3, to_points, 1.0); if (i < (num_fp-4)) { fp += 3; i += 3; } else { break; } } else if (i == (num_fp-3)) { ZnGetBezierPoints(fp, fp+1, fp+1, fp+2, to_points, 1.0); break; } else if (i == (num_fp-2)) { ZnListAdd(to_points, fp+1, ZnListTail); break; } } } /* ********************************************************************************** * * ZnGetCirclePoints -- * Return a pointer to an array of points describing a * circle arc of radius 1.0. The arc is described by start_angle, * end_angle and the type: 0 for arc, 1 for chord, 2 for pie slice, * 3 for a full circle (in which case start_angle and end_angle are * not used. * The number of points is returned in num_points. If type is not 3, * point_list must not be NULL. If not NULL, it is filled with the * computed points. * ********************************************************************************** */ ZnPoint * ZnGetCirclePoints(int type, int quality, ZnReal start_angle, ZnReal angle_extent, unsigned int *num_points, ZnList point_list) { static ZnPoint genarc_finest[] = { /* 128 */ {1.0, 0.0}, {0.99879545617, 0.0490676750517}, {0.99518472653, 0.0980171417729}, {0.989176509646, 0.146730476607}, {0.980785279837, 0.195090324861}, {0.970031252314, 0.24298018342}, {0.956940334469, 0.290284681418}, {0.941544063473, 0.336889858172}, {0.923879530291, 0.382683437725}, {0.903989290333, 0.42755509933}, {0.88192126093, 0.471396743221}, {0.857728605899, 0.514102751035}, {0.831469607468, 0.555570240255}, {0.803207525865, 0.595699312064}, {0.773010446922, 0.634393292011}, {0.74095111805, 0.671558962907}, {0.707106772982, 0.707106789391}, {0.671558945713, 0.740951133634}, {0.634393274074, 0.773010461643}, {0.595699293426, 0.803207539688}, {0.555570220961, 0.83146962036}, {0.514102731131, 0.857728617829}, {0.471396722756, 0.881921271869}, {0.427555078353, 0.903989300254}, {0.382683416286, 0.923879539171}, {0.336889836323, 0.94154407129}, {0.290284659212, 0.956940341205}, {0.242980160911, 0.970031257952}, {0.195090302102, 0.980785284364}, {0.146730453653, 0.98917651305}, {0.0980171186795, 0.995184728805}, {0.0490676518746, 0.998795457308}, {-2.32051033331e-08, 1.0}, {-0.0490676982289, 0.998795455031}, {-0.0980171648663, 0.995184724256}, {-0.146730499561, 0.989176506241}, {-0.19509034762, 0.98078527531}, {-0.24298020593, 0.970031246675}, {-0.290284703624, 0.956940327733}, {-0.33688988002, 0.941544055655}, {-0.382683459163, 0.923879521411}, {-0.427555120307, 0.903989280412}, {-0.471396763686, 0.881921249991}, {-0.514102770939, 0.85772859397}, {-0.555570259549, 0.831469594576}, {-0.595699330703, 0.803207512042}, {-0.634393309949, 0.773010432201}, {-0.6715589801, 0.740951102467}, {-0.707106805799, 0.707106756574}, {-0.740951149217, 0.671558928519}, {-0.773010476365, 0.634393256136}, {-0.803207553511, 0.595699274787}, {-0.831469633252, 0.555570201666}, {-0.857728629759, 0.514102711228}, {-0.881921282808, 0.471396702291}, {-0.903989310176, 0.427555057376}, {-0.923879548052, 0.382683394847}, {-0.941544079108, 0.336889814474}, {-0.956940347941, 0.290284637006}, {-0.97003126359, 0.242980138401}, {-0.980785288892, 0.195090279343}, {-0.989176516455, 0.146730430699}, {-0.995184731079, 0.0980170955862}, {-0.998795458447, 0.0490676286974}, {-1.0, -4.64102066663e-08}, {-0.998795453892, -0.049067721406}, {-0.995184721981, -0.0980171879596}, {-0.989176502836, -0.146730522515}, {-0.980785270783, -0.195090370379}, {-0.970031241037, -0.24298022844}, {-0.956940320997, -0.29028472583}, {-0.941544047838, -0.336889901869}, {-0.923879512531, -0.382683480602}, {-0.90398927049, -0.427555141284}, {-0.881921239052, -0.471396784151}, {-0.85772858204, -0.514102790842}, {-0.831469581684, -0.555570278844}, {-0.803207498218, -0.595699349341}, {-0.77301041748, -0.634393327887}, {-0.740951086883, -0.671558997294}, {-0.707106740165, -0.707106822208}, {-0.671558911325, -0.740951164801}, {-0.634393238198, -0.773010491086}, {-0.595699256149, -0.803207567335}, {-0.555570182372, -0.831469646144}, {-0.514102691324, -0.857728641689}, {-0.471396681826, -0.881921293746}, {-0.427555036399, -0.903989320097}, {-0.382683373409, -0.923879556932}, {-0.336889792626, -0.941544086926}, {-0.2902846148, -0.956940354677}, {-0.242980115891, -0.970031269229}, {-0.195090256583, -0.980785293419}, {-0.146730407745, -0.98917651986}, {-0.0980170724928, -0.995184733354}, {-0.0490676055202, -0.998795459585}, {6.96153097774e-08, -1.0}, {0.0490677445832, -0.998795452754}, {0.098017211053, -0.995184719707}, {0.146730545469, -0.989176499431}, {0.195090393139, -0.980785266256}, {0.242980250949, -0.970031235398}, {0.290284748036, -0.956940314261}, {0.336889923717, -0.94154404002}, {0.382683502041, -0.923879503651}, {0.427555162262, -0.903989260569}, {0.471396804617, -0.881921228114}, {0.514102810746, -0.85772857011}, {0.555570298138, -0.831469568792}, {0.59569936798, -0.803207484395}, {0.634393345825, -0.773010402759}, {0.671559014488, -0.740951071299}, {0.707106838616, -0.707106723757}, {0.740951180385, -0.671558894131}, {0.773010505807, -0.63439322026}, {0.803207581158, -0.59569923751}, {0.831469659036, -0.555570163078}, {0.857728653619, -0.51410267142}, {0.881921304685, -0.471396661361}, {0.903989330019, -0.427555015421}, {0.923879565812, -0.38268335197}, {0.941544094743, -0.336889770777}, {0.956940361414, -0.290284592594}, {0.970031274867, -0.242980093382}, {0.980785297946, -0.195090233824}, {0.989176523265, -0.146730384792}, {0.995184735628, -0.0980170493994}, {0.998795460724, -0.0490675823431}, {1.0, 0.0} }; static ZnPoint genarc_finer[] = { /* 64 */ {1.0, 0.0}, {0.99518472653, 0.0980171417729}, {0.980785279837, 0.195090324861}, {0.956940334469, 0.290284681418}, {0.923879530291, 0.382683437725}, {0.88192126093, 0.471396743221}, {0.831469607468, 0.555570240255}, {0.773010446922, 0.634393292011}, {0.707106772982, 0.707106789391}, {0.634393274074, 0.773010461643}, {0.555570220961, 0.83146962036}, {0.471396722756, 0.881921271869}, {0.382683416286, 0.923879539171}, {0.290284659212, 0.956940341205}, {0.195090302102, 0.980785284364}, {0.0980171186795, 0.995184728805}, {-2.32051033331e-08, 1.0}, {-0.0980171648663, 0.995184724256}, {-0.19509034762, 0.98078527531}, {-0.290284703624, 0.956940327733}, {-0.382683459163, 0.923879521411}, {-0.471396763686, 0.881921249991}, {-0.555570259549, 0.831469594576}, {-0.634393309949, 0.773010432201}, {-0.707106805799, 0.707106756574}, {-0.773010476365, 0.634393256136}, {-0.831469633252, 0.555570201666}, {-0.881921282808, 0.471396702291}, {-0.923879548052, 0.382683394847}, {-0.956940347941, 0.290284637006}, {-0.980785288892, 0.195090279343}, {-0.995184731079, 0.0980170955862}, {-1.0, -4.64102066663e-08}, {-0.995184721981, -0.0980171879596}, {-0.980785270783, -0.195090370379}, {-0.956940320997, -0.29028472583}, {-0.923879512531, -0.382683480602}, {-0.881921239052, -0.471396784151}, {-0.831469581684, -0.555570278844}, {-0.77301041748, -0.634393327887}, {-0.707106740165, -0.707106822208}, {-0.634393238198, -0.773010491086}, {-0.555570182372, -0.831469646144}, {-0.471396681826, -0.881921293746}, {-0.382683373409, -0.923879556932}, {-0.2902846148, -0.956940354677}, {-0.195090256583, -0.980785293419}, {-0.0980170724928, -0.995184733354}, {6.96153097774e-08, -1.0}, {0.098017211053, -0.995184719707}, {0.195090393139, -0.980785266256}, {0.290284748036, -0.956940314261}, {0.382683502041, -0.923879503651}, {0.471396804617, -0.881921228114}, {0.555570298138, -0.831469568792}, {0.634393345825, -0.773010402759}, {0.707106838616, -0.707106723757}, {0.773010505807, -0.63439322026}, {0.831469659036, -0.555570163078}, {0.881921304685, -0.471396661361}, {0.923879565812, -0.38268335197}, {0.956940361414, -0.290284592594}, {0.980785297946, -0.195090233824}, {0.995184735628, -0.0980170493994}, {1.0, 0.0} }; static ZnPoint genarc_fine[] = { /* 40 */ {1.0, 0.0}, {0.987688340232, 0.156434467332}, {0.951056514861, 0.309016998789}, {0.891006521028, 0.453990505942}, {0.809016988919, 0.587785259802}, {0.707106772982, 0.707106789391}, {0.587785241028, 0.809017002559}, {0.453990485266, 0.891006531563}, {0.309016976719, 0.951056522032}, {0.156434444413, 0.987688343862}, {-2.32051033331e-08, 1.0}, {-0.156434490252, 0.987688336602}, {-0.309017020858, 0.95105650769}, {-0.453990526618, 0.891006510493}, {-0.587785278575, 0.809016975279}, {-0.707106805799, 0.707106756574}, {-0.809017016198, 0.587785222255}, {-0.891006542098, 0.453990464591}, {-0.951056529203, 0.30901695465}, {-0.987688347492, 0.156434421493}, {-1.0, -4.64102066663e-08}, {-0.987688332972, -0.156434513171}, {-0.951056500519, -0.309017042928}, {-0.891006499958, -0.453990547294}, {-0.80901696164, -0.587785297348}, {-0.707106740165, -0.707106822208}, {-0.587785203482, -0.809017029838}, {-0.453990443915, -0.891006552633}, {-0.309016932581, -0.951056536373}, {-0.156434398574, -0.987688351122}, {6.96153097774e-08, -1.0}, {0.15643453609, -0.987688329342}, {0.309017064997, -0.951056493349}, {0.45399056797, -0.891006489423}, {0.587785316122, -0.809016948}, {0.707106838616, -0.707106723757}, {0.809017043478, -0.587785184709}, {0.891006563167, -0.453990423239}, {0.951056543544, -0.309016910511}, {0.987688354752, -0.156434375655}, {1.0, 0.0} }; static ZnPoint genarc_medium[] = { /* 20 */ {1.0, 0.0}, {0.951056514861, 0.309016998789}, {0.809016988919, 0.587785259802}, {0.587785241028, 0.809017002559}, {0.309016976719, 0.951056522032}, {-2.32051033331e-08, 1.0}, {-0.309017020858, 0.95105650769}, {-0.587785278575, 0.809016975279}, {-0.809017016198, 0.587785222255}, {-0.951056529203, 0.30901695465}, {-1.0, -4.64102066663e-08}, {-0.951056500519, -0.309017042928}, {-0.80901696164, -0.587785297348}, {-0.587785203482, -0.809017029838}, {-0.309016932581, -0.951056536373}, {6.96153097774e-08, -1.0}, {0.309017064997, -0.951056493349}, {0.587785316122, -0.809016948}, {0.809017043478, -0.587785184709}, {0.951056543544, -0.309016910511}, {1.0, 0.0} }; static ZnPoint genarc_coarse[] = { /* 10 */ {1.0, 0.0}, {0.809016988919, 0.587785259802}, {0.309016976719, 0.951056522032}, {-0.309017020858, 0.95105650769}, {-0.809017016198, 0.587785222255}, {-1.0, -4.64102066663e-08}, {-0.80901696164, -0.587785297348}, {-0.309016932581, -0.951056536373}, {0.309017064997, -0.951056493349}, {0.809017043478, -0.587785184709}, {1.0, 0.0} }; unsigned int num_p, i; ZnPoint *p, *p_from; ZnPoint center_p = { 0.0, 0.0 }; ZnPoint start_p, wp; ZnReal iangle, end_angle=0; switch (quality) { case ZN_CIRCLE_COARSE: num_p = sizeof(genarc_coarse)/sizeof(ZnPoint); p = p_from = genarc_coarse; break; case ZN_CIRCLE_MEDIUM: num_p = sizeof(genarc_medium)/sizeof(ZnPoint); p = p_from = genarc_medium; break; case ZN_CIRCLE_FINER: num_p = sizeof(genarc_finer)/sizeof(ZnPoint); p = p_from = genarc_finer; break; case ZN_CIRCLE_FINEST: num_p = sizeof(genarc_finest)/sizeof(ZnPoint); p = p_from = genarc_finest; break; default: case ZN_CIRCLE_FINE: num_p = sizeof(genarc_fine)/sizeof(ZnPoint); p = p_from = genarc_fine; } /*printf("start: %g, extent: %g\n", start_angle, angle_extent);*/ if (angle_extent == 2*M_PI) { type = 3; } if (type != 3) { end_angle = start_angle+angle_extent; if (angle_extent < 0) { iangle = start_angle; start_angle = end_angle; end_angle = iangle; } /* * normalize start_angle and end_angle. */ if (start_angle < 0.0) { start_angle += 2.0*M_PI; } if (end_angle < 0.0) { end_angle += 2.0*M_PI; } if (end_angle < start_angle) { end_angle += 2.0*M_PI; } /*printf("---start: %g, end: %g\n", start_angle, end_angle);*/ } /* * Now 0 <= start_angle < 2 * M_PI and start_angle <= end_angle. */ if ((type != 3) || (point_list != NULL)) { if (type == 3) { ZnListAssertSize(point_list, num_p); p = ZnListArray(point_list); for (i = 0; i < num_p; i++, p++, p_from++) { *p = *p_from; } } else { ZnListEmpty(point_list); iangle = 2*M_PI / (num_p-1); start_p.x = cos(start_angle); start_p.y = sin(start_angle); ZnListAdd(point_list, &start_p, ZnListTail); i = (unsigned int) (start_angle / iangle); if ((i * iangle) < start_angle) { i++; } while (1) { if (start_angle + iangle <= end_angle) { if (i == num_p-1) { i = 0; } ZnListAdd(point_list, &p_from[i], ZnListTail); start_angle += iangle; i++; } else { wp.x = cos(end_angle); wp.y = sin(end_angle); ZnListAdd(point_list, &wp, ZnListTail); break; } } if (type == 1) { ZnListAdd(point_list, &start_p, ZnListTail); } else if (type == 2) { ZnListAdd(point_list, ¢er_p, ZnListTail); ZnListAdd(point_list, &start_p, ZnListTail); } } p = ZnListArray(point_list); num_p = ZnListSize(point_list); } *num_points = num_p; return p; } /* ********************************************************************************** * * ZnGetArcPath -- * Compute in to_points a set of Bezier control points describing an arc * path given the start angle, the stop angle and the type: 0 for arc, * 1 for chord, 2 for pie slice. * To obtain the actual polygonal shape, the client should use GetBezierPath * on the returned controls (after applying transform). The returned arc * is circular and centered on 0,0. * ********************************************************************************** */ static ZnReal arc_nodes_x[4] = { 1.0, 0.0, -1.0, 0.0 }; static ZnReal arc_nodes_y[4] = { 0.0, 1.0, 0.0, -1.0 }; static ZnReal arc_controls_x[8] = { 1.0, 0.55197, -0.55197, -1.0, -1.0, -0.55197, 0.55197, 1.0 }; static ZnReal arc_controls_y[8] = { 0.55197, 1.0, 1.0, 0.55197, -0.55197, -1.0, -1.0, -0.55197 }; void ZnGetArcPath(ZnReal start_angle, ZnReal end_angle, int type, ZnList to_points) { int start_quad, end_quad, quadrant; ZnPoint center_p = { 0.0, 0.0 }; ZnPoint start_p = center_p; /* * make sure the output vector is empty. */ ZnListEmpty(to_points); /* * normalize start_angle and end_angle. */ start_angle = fmod(start_angle, 2.0*M_PI); if (start_angle < 0.0) { start_angle += 2.0*M_PI; } end_angle = fmod(end_angle, 2.0*M_PI); if (end_angle < 0.0) { end_angle += 2.0*M_PI; } if (start_angle >= end_angle) { if (start_angle == end_angle) { type = 3; } end_angle += 2.0*M_PI; } /* * Now 0 <= start_angle < 2 * M_PI and start_angle <= end_angle. */ start_quad = (int) (start_angle / (M_PI / 2.0)); end_quad = (int) (end_angle / (M_PI / 2.0)); for (quadrant = start_quad; quadrant <= end_quad; quadrant++) { ZnPoint controls[4]; ZnReal t; controls[0].x = arc_nodes_x[quadrant % 4]; controls[0].y = arc_nodes_y[quadrant % 4]; controls[1].x = arc_controls_x[2 * (quadrant % 4)]; controls[1].y = arc_controls_y[2 * (quadrant % 4)]; controls[2].x = arc_controls_x[2 * (quadrant % 4) + 1]; controls[2].y = arc_controls_y[2 * (quadrant % 4) + 1]; controls[3].x = arc_nodes_x[(quadrant + 1) % 4]; controls[3].y = arc_nodes_y[(quadrant + 1) % 4]; if (quadrant == start_quad) { t = Arc2Param(controls, start_angle); BezierSubdivide(controls, t, False); /* * The path is still empty and we have to create the first * vertex. */ start_p = controls[0]; ZnListAdd(to_points, &controls[0], ZnListTail); } if (quadrant == end_quad) { t = Arc2Param(controls, end_angle); if (!t) { break; } BezierSubdivide(controls, t, True); } ZnListAdd(to_points, &controls[1], ZnListTail); ZnListAdd(to_points, &controls[2], ZnListTail); ZnListAdd(to_points, &controls[3], ZnListTail); } if (type == 2) { ZnListAdd(to_points, ¢er_p, ZnListTail); /* * Can't add this one, it would be interpreted by GetBezierPath * as an off-curve control. The path should be closed by the client * after expansion by GetBezierPath. * ZnListAdd(to_points, &start_p, ZnListTail); */ } else if (type == 1) { ZnListAdd(to_points, &start_p, ZnListTail); } } /* ********************************************************************************** * * SmoothPathWithBezier -- * Compute in to_points a new set of points describing a smoothed path based * on the path given in from_points. The algorithm use Bezier cubic curves. * ********************************************************************************** */ void ZnSmoothPathWithBezier(ZnPoint *fp, unsigned int num_fp, ZnList to_points) { ZnBool closed; ZnPoint s[4]; unsigned int i; /* * make sure the output vector is empty */ ZnListEmpty(to_points); /* * If the curve is closed, generates a Bezier curve that * spans the closure. Else simply add the first point to * the path. */ if ((fp[0].x == fp[num_fp-1].x) && (fp[0].y == fp[num_fp-1].y)) { closed = True; s[0].x = 0.5*fp[num_fp-2].x + 0.5*fp[0].x; s[0].y = 0.5*fp[num_fp-2].y + 0.5*fp[0].y; s[1].x = 0.167*fp[num_fp-2].x + 0.833*fp[0].x; s[1].y = 0.167*fp[num_fp-2].y + 0.833*fp[0].y; s[2].x = 0.833*fp[0].x + 0.167*fp[1].x; s[2].y = 0.833*fp[0].y + 0.167*fp[1].y; s[3].x = 0.5*fp[0].x + 0.5*fp[1].x; s[3].y = 0.5*fp[0].y + 0.5*fp[1].y; ZnListAdd(to_points, s, ZnListTail); ZnGetBezierPoints(s, s+1, s+2, s+3, to_points, 1.0); } else { closed = False; ZnListAdd(to_points, &fp[0], ZnListTail); } for (i = 2; i < num_fp; i++, fp++) { /* * Setup the first two control points. This differ * for first segment of open curves. */ if ((i == 2) && !closed) { s[0] = fp[0]; s[1].x = 0.333*fp[0].x + 0.667*fp[1].x; s[1].y = 0.333*fp[0].y + 0.667*fp[1].y; } else { s[0].x = 0.5*fp[0].x + 0.5*fp[1].x; s[0].y = 0.5*fp[0].y + 0.5*fp[1].y; s[1].x = 0.167*fp[0].x + 0.833*fp[1].x; s[1].y = 0.167*fp[0].y + 0.833*fp[1].y; } /* * Setup the last two control points. This also differ * for last segment of open curves. */ if ((i == num_fp-1) && !closed) { s[2].x = 0.667*fp[1].x + 0.333*fp[2].x; s[2].y = 0.667*fp[1].y + 0.333*fp[2].y; s[3] = fp[2]; } else { s[2].x = 0.833*fp[1].x + 0.167*fp[2].x; s[2].y = 0.833*fp[1].y + 0.167*fp[2].y; s[3].x = 0.5*fp[1].x + 0.5*fp[2].x; s[3].y = 0.5*fp[1].y + 0.5*fp[2].y; } /* * If the first two points or the last two are equal * output the last control point. Else generate the * Bezier curve. */ if (((fp[0].x == fp[1].x) && (fp[0].y == fp[1].y)) || ((fp[1].x == fp[2].x) && (fp[1].y == fp[2].y))) { ZnListAdd(to_points, &s[3], ZnListTail); } else { ZnGetBezierPoints(s, s+1, s+2, s+3, to_points, 1.0); } } } /* ********************************************************************************** * * FitBezier -- * Fit a Bezier curve to a (sub)set of digitized points. * * From: An Algorithm for Automatically Fitting Digitized Curves * by Philip J. Schneider in "Graphics Gems", Academic Press, 1990 * ********************************************************************************** */ static ZnReal V2DistanceBetween2Points(ZnPoint *a, ZnPoint *b) { ZnReal dx = a->x - b->x; ZnReal dy = a->y - b->y; return sqrt((dx*dx)+(dy*dy)); } static ZnReal V2SquaredLength(ZnPoint *a) { return (a->x * a->x)+(a->y * a->y); } static ZnReal V2Length(ZnPoint *a) { return sqrt(V2SquaredLength(a)); } static ZnPoint * V2Scale(ZnPoint *v, ZnReal newlen) { ZnReal len = V2Length(v); if (len != 0.0) { v->x *= newlen/len; v->y *= newlen/len; } return v; } static ZnPoint * V2Negate(ZnPoint *v) { v->x = -v->x; v->y = -v->y; return v; } static ZnPoint * V2Normalize(ZnPoint *v) { ZnReal len = sqrt(V2Length(v)); if (len != 0.0) { v->x /= len; v->y /= len; } return v; } static ZnPoint * V2Add(ZnPoint *a, ZnPoint *b, ZnPoint *c) { c->x = a->x + b->x; c->y = a->y + b->y; return c; } static ZnReal V2Dot(ZnPoint *a, ZnPoint *b) { return (a->x*b->x) + (a->y*b->y); } static ZnPoint V2AddII(ZnPoint a, ZnPoint b) { ZnPoint c; c.x = a.x + b.x; c.y = a.y + b.y; return c; } static ZnPoint V2ScaleIII(ZnPoint v, ZnReal s) { ZnPoint result; result.x = v.x * s; result.y = v.y * s; return result; } static ZnPoint V2SubII(ZnPoint a, ZnPoint b) { ZnPoint c; c.x = a.x - b.x; c.y = a.y - b.y; return c; } /* * B0, B1, B2, B3, Bezier multipliers. */ static ZnReal B0(ZnReal u) { ZnReal tmp = 1.0 - u; return tmp * tmp * tmp; } static ZnReal B1(ZnReal u) { ZnReal tmp = 1.0 - u; return 3 * u * (tmp * tmp); } static ZnReal B2(ZnReal u) { ZnReal tmp = 1.0 - u; return 3 * u * u * tmp; } static ZnReal B3(ZnReal u) { return u * u * u; } /* * ChordLengthParameterize -- * Assign parameter values to digitized points * using relative distances between points. */ static ZnReal * ChordLengthParameterize(ZnPoint *d, unsigned int first, unsigned int last) { unsigned int i; ZnReal *u; u = (ZnReal *) ZnMalloc((unsigned) (last-first+1) * sizeof(ZnReal)); u[0] = 0.0; for (i = first+1; i <= last; i++) { u[i-first] = u[i-first-1] + V2DistanceBetween2Points(&d[i], &d[i-1]); } for (i = first + 1; i <= last; i++) { u[i-first] = u[i-first] / u[last-first]; } return u; } /* * Bezier -- * Evaluate a Bezier curve at a particular parameter value * */ static ZnPoint BezierII(int degree, ZnPoint *V, ZnReal t) { int i, j; ZnPoint Q; /* Point on curve at parameter t */ ZnPoint *Vtemp; /* Local copy of control points */ /* Copy array */ Vtemp = (ZnPoint *) ZnMalloc((unsigned)((degree+1) * sizeof (ZnPoint))); for (i = 0; i <= degree; i++) { Vtemp[i] = V[i]; } /* Triangle computation */ for (i = 1; i <= degree; i++) { for (j = 0; j <= degree-i; j++) { Vtemp[j].x = (1.0 - t) * Vtemp[j].x + t * Vtemp[j+1].x; Vtemp[j].y = (1.0 - t) * Vtemp[j].y + t * Vtemp[j+1].y; } } Q = Vtemp[0]; ZnFree(Vtemp); return Q; } /* * NewtonRaphsonRootFind -- * Use Newton-Raphson iteration to find better root. */ static ZnReal NewtonRaphsonRootFind(ZnPoint *Q, ZnPoint P, ZnReal u) { ZnReal numerator, denominator; ZnPoint Q1[3], Q2[2]; /* Q' and Q'' */ ZnPoint Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */ ZnReal uPrime; /* Improved u */ unsigned int i; /* Compute Q(u) */ Q_u = BezierII(3, Q, u); /* Generate control vertices for Q' */ for (i = 0; i <= 2; i++) { Q1[i].x = (Q[i+1].x - Q[i].x) * 3.0; Q1[i].y = (Q[i+1].y - Q[i].y) * 3.0; } /* Generate control vertices for Q'' */ for (i = 0; i <= 1; i++) { Q2[i].x = (Q1[i+1].x - Q1[i].x) * 2.0; Q2[i].y = (Q1[i+1].y - Q1[i].y) * 2.0; } /* Compute Q'(u) and Q''(u) */ Q1_u = BezierII(2, Q1, u); Q2_u = BezierII(1, Q2, u); /* Compute f(u)/f'(u) */ numerator = (Q_u.x - P.x) * (Q1_u.x) + (Q_u.y - P.y) * (Q1_u.y); denominator = (Q1_u.x) * (Q1_u.x) + (Q1_u.y) * (Q1_u.y) + (Q_u.x - P.x) * (Q2_u.x) + (Q_u.y - P.y) * (Q2_u.y); /* u = u - f(u)/f'(u) */ uPrime = u - (numerator/denominator); return (uPrime); } /* * Reparameterize -- * Given set of points and their parameterization, try to find * a better parameterization. */ static ZnReal * Reparameterize(ZnPoint *d, unsigned int first, unsigned int last, ZnReal *u, ZnPoint *bezCurve) { unsigned int nPts = last-first+1; unsigned int i; ZnReal *uPrime; /* New parameter values */ uPrime = (ZnReal *) ZnMalloc(nPts * sizeof(ZnReal)); for (i = first; i <= last; i++) { uPrime[i-first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i-first]); } return (uPrime); } /* * GenerateBezier -- * Use least-squares method to find Bezier control * points for region. */ static void GenerateBezier(ZnPoint *d, unsigned int first, unsigned int last, ZnReal *uPrime, ZnPoint tHat1, ZnPoint tHat2, ZnPoint *bez_curve) { unsigned int i; ZnPoint *A0, *A1; /* Precomputed rhs for eqn */ unsigned int num_points; /* Number of pts in sub-curve */ ZnReal C[2][2]; /* Matrix C */ ZnReal X[2]; /* Matrix X */ ZnReal det_C0_C1; /* Determinants of matrices */ ZnReal det_C0_X, det_X_C1; ZnReal alpha_l; /* Alpha values, left and right */ ZnReal alpha_r; ZnPoint tmp; /* Utility variable */ num_points = last - first + 1; A0 = (ZnPoint *) ZnMalloc(num_points * sizeof(ZnPoint)); A1 = (ZnPoint *) ZnMalloc(num_points * sizeof(ZnPoint)); /* Compute the A's */ for (i = 0; i < num_points; i++) { ZnPoint v1, v2; v1 = tHat1; v2 = tHat2; V2Scale(&v1, B1(uPrime[i])); V2Scale(&v2, B2(uPrime[i])); A0[i] = v1; A1[i] = v2; } /* Create the C and X matrices */ C[0][0] = 0.0; C[0][1] = 0.0; C[1][0] = 0.0; C[1][1] = 0.0; X[0] = 0.0; X[1] = 0.0; for (i = 0; i < num_points; i++) { C[0][0] += V2Dot(&A0[i], &A0[i]); C[0][1] += V2Dot(&A0[i], &A1[i]); C[1][0] = C[0][1]; C[1][1] += V2Dot(&A1[i], &A1[i]); tmp = V2SubII(d[first + i], V2AddII(V2ScaleIII(d[first], B0(uPrime[i])), V2AddII(V2ScaleIII(d[first], B1(uPrime[i])), V2AddII(V2ScaleIII(d[last], B2(uPrime[i])), V2ScaleIII(d[last], B3(uPrime[i])))))); X[0] += V2Dot(&A0[i], &tmp); X[1] += V2Dot(&A1[i], &tmp); } /* Compute the determinants of C and X */ det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1]; det_C0_X = C[0][0] * X[1] - C[0][1] * X[0]; det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1]; /* Finally, derive alpha values */ if (det_C0_C1 == 0.0) { det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12; } alpha_l = det_X_C1 / det_C0_C1; alpha_r = det_C0_X / det_C0_C1; /* If alpha negative, use the Wu/Barsky heuristic (see text) */ if (alpha_l < 0.0 || alpha_r < 0.0) { ZnReal dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0; bez_curve[0] = d[first]; bez_curve[3] = d[last]; V2Add(&bez_curve[0], V2Scale(&tHat1, dist), &bez_curve[1]); V2Add(&bez_curve[3], V2Scale(&tHat2, dist), &bez_curve[2]); } else { /* First and last control points of the Bezier curve are */ /* positioned exactly at the first and last data points */ /* Control points 1 and 2 are positioned an alpha distance out */ /* on the tangent vectors, left and right, respectively */ bez_curve[0] = d[first]; bez_curve[3] = d[last]; V2Add(&bez_curve[0], V2Scale(&tHat1, alpha_l), &bez_curve[1]); V2Add(&bez_curve[3], V2Scale(&tHat2, alpha_r), &bez_curve[2]); } ZnFree(A0); ZnFree(A1); } /* * ComputeMaxError -- * Find the maximum squared distance of digitized points * to fitted curve. */ static ZnReal ComputeMaxError(ZnPoint *d, unsigned int first, unsigned int last, ZnPoint *bez_curve, ZnReal *u, unsigned int *splitPoint) { unsigned int i; ZnReal maxDist; /* Maximum error */ ZnReal dist; /* Current error */ ZnPoint P; /* Point on curve */ ZnPoint v; /* Vector from point to curve */ *splitPoint = (last - first + 1)/2; maxDist = 0.0; for (i = first + 1; i < last; i++) { P = BezierII(3, bez_curve, u[i-first]); v = V2SubII(P, d[i]); dist = V2SquaredLength(&v); if (dist >= maxDist) { maxDist = dist; *splitPoint = i; } } return (maxDist); } /* * ComputeLeftTangent, * ComputeRightTangent, * ComputeCenterTangent -- * Approximate unit tangents at endpoints and * center of digitized curve. */ static ZnPoint ComputeLeftTangent(ZnPoint *d, unsigned int end) { ZnPoint tHat1; tHat1 = V2SubII(d[end+1], d[end]); tHat1 = *V2Normalize(&tHat1); return tHat1; } static ZnPoint ComputeRightTangent(ZnPoint *d, unsigned int end) { ZnPoint tHat2; tHat2 = V2SubII(d[end-1], d[end]); tHat2 = *V2Normalize(&tHat2); return tHat2; } static ZnPoint ComputeCenterTangent(ZnPoint *d, unsigned int center) { ZnPoint V1, V2, tHatCenter; V1 = V2SubII(d[center-1], d[center]); V2 = V2SubII(d[center], d[center+1]); tHatCenter.x = (V1.x + V2.x)/2.0; tHatCenter.y = (V1.y + V2.y)/2.0; tHatCenter = *V2Normalize(&tHatCenter); return tHatCenter; } static void FitCubic(ZnPoint *d, unsigned int first, unsigned int last, ZnPoint tHat1, ZnPoint tHat2, ZnReal error, ZnList controls) { ZnPoint *bez_curve; /* Control points of fitted Bezier curve*/ ZnReal *u; /* Parameter values for point */ ZnReal *uPrime; /* Improved parameter values */ ZnReal max_err; /* Maximum fitting error */ unsigned int splitPoint; /* Point to split point set at */ unsigned int num_points; /* Number of points in subset */ ZnReal iteration_err; /* Error below which you try iterating */ unsigned int max_iter = 4; /* Max times to try iterating */ ZnPoint tHatCenter; /* Unit tangent vector at splitPoint */ unsigned int i; iteration_err = error * error; num_points = last - first + 1; ZnListAssertSize(controls, ZnListSize(controls)+4); bez_curve = ZnListAt(controls, ZnListSize(controls)-4); /* Use heuristic if region only has two points in it */ if (num_points == 2) { ZnReal dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0; bez_curve[0] = d[first]; bez_curve[3] = d[last]; V2Add(&bez_curve[0], V2Scale(&tHat1, dist), &bez_curve[1]); V2Add(&bez_curve[3], V2Scale(&tHat2, dist), &bez_curve[2]); return; } /* Parameterize points, and attempt to fit curve */ u = ChordLengthParameterize(d, first, last); GenerateBezier(d, first, last, u, tHat1, tHat2, bez_curve); /* Find max deviation of points to fitted curve */ max_err = ComputeMaxError(d, first, last, bez_curve, u, &splitPoint); if (max_err < error) { ZnFree(u); return; } /* * If error not too large, try some reparameterization * and iteration. */ if (max_err < iteration_err) { for (i = 0; i < max_iter; i++) { uPrime = Reparameterize(d, first, last, u, bez_curve); GenerateBezier(d, first, last, uPrime, tHat1, tHat2, bez_curve); max_err = ComputeMaxError(d, first, last, bez_curve, uPrime, &splitPoint); if (max_err < error) { ZnFree(u); return; } ZnFree(u); u = uPrime; } } /* Fitting failed -- split at max error point and fit recursively */ ZnFree(u); ZnListAssertSize(controls, ZnListSize(controls)-4); tHatCenter = ComputeCenterTangent(d, splitPoint); FitCubic(d, first, splitPoint, tHat1, tHatCenter, error, controls); V2Negate(&tHatCenter); FitCubic(d, splitPoint, last, tHatCenter, tHat2, error, controls); } void ZnFitBezier(ZnPoint *pts, unsigned int num_points, ZnReal error, ZnList controls) { ZnPoint tHat1, tHat2; /* Unit tangent vectors at endpoints */ tHat1 = ComputeLeftTangent(pts, 0); tHat2 = ComputeRightTangent(pts, num_points-1); FitCubic(pts, 0, num_points-1, tHat1, tHat2, error, controls); }