From 4c9faffa908a0749968aed6638f6f54d42e1d3cb Mon Sep 17 00:00:00 2001 From: lecoanet Date: Thu, 11 May 2000 14:05:15 +0000 Subject: Ajout des fonctions de manipulation de la structure de polygone. Exportation de la fonction de lissage par des Beziers. Ajout de la fonction d'interpolation par des Beziers. --- generic/Geo.c | 618 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 608 insertions(+), 10 deletions(-) (limited to 'generic/Geo.c') diff --git a/generic/Geo.c b/generic/Geo.c index c564be7..aca1590 100644 --- a/generic/Geo.c +++ b/generic/Geo.c @@ -35,12 +35,76 @@ #include "WidgetInfo.h" #include +#include static const char rcsid[] = "$Id$"; static const char compile_id[]="$Compile: " __FILE__ " " __DATE__ " " __TIME__ " $"; +void +POLY_INIT(ZnPoly *poly) { + poly->num_contours = 0; + poly->contours = NULL; + poly->holes = NULL; + poly->cw = NULL; +} + +void +POLY_CONTOUR1(ZnPoly *poly, + ZnPoint *pts, + int num_pts) { + poly->num_contours = 1; + poly->holes = &poly->hole1; + poly->hole1 = False; + poly->cw = &poly->cw1; + poly->contours = &poly->contour1; + poly->contour1.num_points = num_pts; + poly->contour1.points = pts; +} + +void +POLY_SET(ZnPoly *poly1, + ZnPoly *poly2) { + POLY_FREE(poly1); + if (poly2->num_contours == 1) { + POLY_CONTOUR1(poly1, poly2->contours[0].points, poly2->contours[0].num_points); + if (poly2->contours != &poly2->contour1) { + ZnFree(poly2->contours); + } + } + else { + poly1->num_contours = poly2->num_contours; + poly1->contours = poly2->contours; + poly1->holes = poly2->holes; + poly1->cw = poly2->cw; + } +} + +void +POLY_FREE(ZnPoly *poly) { + if (poly->num_contours) { + int i; + for (i = 0; i < poly->num_contours; i++) { + ZnFree(poly->contours[i].points); + } + if ((poly->contours != &poly->contour1) && + (poly->num_contours > 1)) { + ZnFree(poly->contours); + if (poly->holes) { + ZnFree(poly->holes); + } + if (poly->cw) { + ZnFree(poly->cw); + } + } + poly->num_contours = 0; + poly->contours = NULL; + poly->holes = NULL; + poly->cw = NULL; + } +} + /* * Compute the origin of the rectangle given * by position, anchor, width and height. @@ -828,7 +892,7 @@ PolylineInBBox(ZnPoint *points, * poly[2]. */ if ((join_style == JoinBevel) || do_miter_as_bevel) { - if (PolygonInBBox(poly, 4, bbox) != inside) { + if (PolygonInBBox(poly, 4, bbox, NULL) != inside) { return 0; } do_miter_as_bevel = False; @@ -853,7 +917,7 @@ PolylineInBBox(ZnPoint *points, GetButtPoints(points, &points[1], width, 0, &poly[2], &poly[3]); } - if (PolygonInBBox(poly, 4, bbox) != inside) { + if (PolygonInBBox(poly, 4, bbox, NULL) != inside) { return 0; } } @@ -874,17 +938,23 @@ PolylineInBBox(ZnPoint *points, /* * 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 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 PolygonInBBox(ZnPoint *points, int num_points, - ZnBBox *bbox) + 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; @@ -937,6 +1007,9 @@ PolygonInBBox(ZnPoint *points, num_points, bbox->orig.x, bbox->orig.y, PolygonToPointDist(points, num_points, &bbox->orig));*/ if (PolygonToPointDist(points, num_points, &bbox->orig) <= 0.0) { + if (area_enclosed) { + *area_enclosed = True; + } return 0; } @@ -2058,18 +2131,14 @@ GetArcPath(ZnReal start_angle, ********************************************************************************** */ void -SmoothPathWithBezier(ZnList from_points, +SmoothPathWithBezier(ZnPoint *fp, + int num_fp, ZnList to_points) { - ZnPoint *fp; - int num_fp; ZnBool closed; ZnPoint s[4]; int i; - fp = (ZnPoint *) ZnListArray(from_points); - num_fp = ZnListSize(from_points); - /* * make sure the output vector is empty */ @@ -2150,6 +2219,535 @@ SmoothPathWithBezier(ZnList from_points, /* ********************************************************************************** * + * 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, + int first, + int last) +{ + 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 */ + 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, + int first, + int last, + ZnReal *u, + ZnPoint *bezCurve) +{ + int nPts = last-first+1; + 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, + int first, + int last, + ZnReal *uPrime, + ZnPoint tHat1, + ZnPoint tHat2, + ZnPoint *bez_curve) +{ + int i; + ZnPoint *A0, *A1; /* Precomputed rhs for eqn */ + 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, + int first, + int last, + ZnPoint *bez_curve, + ZnReal *u, + int *splitPoint) +{ + 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, + int end) +{ + ZnPoint tHat1; + tHat1 = V2SubII(d[end+1], d[end]); + tHat1 = *V2Normalize(&tHat1); + return tHat1; +} + +static ZnPoint +ComputeRightTangent(ZnPoint *d, + int end) +{ + ZnPoint tHat2; + tHat2 = V2SubII(d[end-1], d[end]); + tHat2 = *V2Normalize(&tHat2); + return tHat2; +} + + +static ZnPoint +ComputeCenterTangent(ZnPoint *d, + 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, + int first, + 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 */ + int splitPoint; /* Point to split point set at */ + int num_points; /* Number of points in subset */ + ZnReal iteration_err; /* Error below which you try iterating */ + int max_iter = 4; /* Max times to try iterating */ + ZnPoint tHatCenter; /* Unit tangent vector at splitPoint */ + 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 +FitBezier(ZnPoint *pts, + 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); +} + + +/* + ********************************************************************************** + * * GetLineEnd -- * Compute the points describing the given line end style at point p1 for * the line p1,p2. Point p1 is adjusted to fit the line end. -- cgit v1.1