/* Solving the Nearest Point-on-Curve Problem and A Bezier Curve-Based Root-Finder by Philip J. Schneider from "Graphics Gems", Academic Press, 1990 */ /* point_on_curve.c */ #include #include #include #include "GraphicsGems.h" #define TESTMODE /* * Forward declarations */ Point2 NearestPointOnCurve(); static int FindRoots(); static Point2 *ConvertToBezierForm(); static double ComputeXIntercept(); static int ControlPolygonFlatEnough(); static int CrossingCount(); static Point2 Bezier(); static Vector2 V2ScaleII(); int MAXDEPTH = 64; /* Maximum depth for recursion */ #define EPSILON (ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */ #define DEGREE 3 /* Cubic Bezier curve */ #define W_DEGREE 5 /* Degree of eqn to find roots of */ #ifdef TESTMODE /* * main : * Given a cubic Bezier curve (i.e., its control points), and some * arbitrary point in the plane, find the point on the curve * closest to that arbitrary point. */ main() { static Point2 bezCurve[4] = { /* A cubic Bezier curve */ { 0.0, 0.0 }, { 1.0, 2.0 }, { 3.0, 3.0 }, { 4.0, 2.0 }, }; static Point2 arbPoint = { 3.5, 2.0 }; /*Some arbitrary point*/ Point2 pointOnCurve; /* Nearest point on the curve */ /* Find the closest point */ pointOnCurve = NearestPointOnCurve(arbPoint, bezCurve); printf("pointOnCurve : (%4.4f, %4.4f)\n", pointOnCurve.x, pointOnCurve.y); } #endif /* TESTMODE */ /* * NearestPointOnCurve : * Compute the parameter value of the point on a Bezier * curve segment closest to some arbtitrary, user-input point. * Return the point on the curve at that parameter value. * */ Point2 NearestPointOnCurve(P, V) Point2 P; /* The user-supplied point */ Point2 *V; /* Control points of cubic Bezier */ { Point2 *w; /* Ctl pts for 5th-degree eqn */ double t_candidate[W_DEGREE]; /* Possible roots */ int n_solutions; /* Number of roots found */ double t; /* Parameter value of closest pt*/ /* Convert problem to 5th-degree Bezier form */ w = ConvertToBezierForm(P, V); /* Find all possible roots of 5th-degree equation */ n_solutions = FindRoots(w, W_DEGREE, t_candidate, 0); free((char *)w); /* Compare distances of P to all candidates, and to t=0, and t=1 */ { double dist, new_dist; Point2 p; Vector2 v; int i; /* Check distance to beginning of curve, where t = 0 */ dist = V2SquaredLength(V2Sub(&P, &V[0], &v)); t = 0.0; /* Find distances for candidate points */ for (i = 0; i < n_solutions; i++) { p = Bezier(V, DEGREE, t_candidate[i], (Point2 *)NULL, (Point2 *)NULL); new_dist = V2SquaredLength(V2Sub(&P, &p, &v)); if (new_dist < dist) { dist = new_dist; t = t_candidate[i]; } } /* Finally, look at distance to end point, where t = 1.0 */ new_dist = V2SquaredLength(V2Sub(&P, &V[DEGREE], &v)); if (new_dist < dist) { dist = new_dist; t = 1.0; } } /* Return the point on the curve at parameter value t */ printf("t : %4.12f\n", t); return (Bezier(V, DEGREE, t, (Point2 *)NULL, (Point2 *)NULL)); } /* * ConvertToBezierForm : * Given a point and a Bezier curve, generate a 5th-degree * Bezier-format equation whose solution finds the point on the * curve nearest the user-defined point. */ static Point2 *ConvertToBezierForm(P, V) Point2 P; /* The point to find t for */ Point2 *V; /* The control points */ { int i, j, k, m, n, ub, lb; int row, column; /* Table indices */ Vector2 c[DEGREE+1]; /* V(i)'s - P */ Vector2 d[DEGREE]; /* V(i+1) - V(i) */ Point2 *w; /* Ctl pts of 5th-degree curve */ double cdTable[3][4]; /* Dot product of c, d */ static double z[3][4] = { /* Precomputed "z" for cubics */ {1.0, 0.6, 0.3, 0.1}, {0.4, 0.6, 0.6, 0.4}, {0.1, 0.3, 0.6, 1.0}, }; /*Determine the c's -- these are vectors created by subtracting*/ /* point P from each of the control points */ for (i = 0; i <= DEGREE; i++) { V2Sub(&V[i], &P, &c[i]); } /* Determine the d's -- these are vectors created by subtracting*/ /* each control point from the next */ for (i = 0; i <= DEGREE - 1; i++) { d[i] = V2ScaleII(V2Sub(&V[i+1], &V[i], &d[i]), 3.0); } /* Create the c,d table -- this is a table of dot products of the */ /* c's and d's */ for (row = 0; row <= DEGREE - 1; row++) { for (column = 0; column <= DEGREE; column++) { cdTable[row][column] = V2Dot(&d[row], &c[column]); } } /* Now, apply the z's to the dot products, on the skew diagonal*/ /* Also, set up the x-values, making these "points" */ w = (Point2 *)malloc((unsigned)(W_DEGREE+1) * sizeof(Point2)); for (i = 0; i <= W_DEGREE; i++) { w[i].y = 0.0; w[i].x = (double)(i) / W_DEGREE; } n = DEGREE; m = DEGREE-1; for (k = 0; k <= n + m; k++) { lb = MAX(0, k - m); ub = MIN(k, n); for (i = lb; i <= ub; i++) { j = k - i; w[i+j].y += cdTable[j][i] * z[j][i]; } } return (w); } /* * FindRoots : * Given a 5th-degree equation in Bernstein-Bezier form, find * all of the roots in the interval [0, 1]. Return the number * of roots found. */ static int FindRoots(w, degree, t, depth) Point2 *w; /* The control points */ int degree; /* The degree of the polynomial */ double *t; /* RETURN candidate t-values */ int depth; /* The depth of the recursion */ { int i; Point2 Left[W_DEGREE+1], /* New left and right */ Right[W_DEGREE+1]; /* control polygons */ int left_count, /* Solution count from */ right_count; /* children */ double left_t[W_DEGREE+1], /* Solutions from kids */ right_t[W_DEGREE+1]; switch (CrossingCount(w, degree)) { case 0 : { /* No solutions here */ return 0; } case 1 : { /* Unique solution */ /* Stop recursion when the tree is deep enough */ /* if deep enough, return 1 solution at midpoint */ if (depth >= MAXDEPTH) { t[0] = (w[0].x + w[W_DEGREE].x) / 2.0; return 1; } if (ControlPolygonFlatEnough(w, degree)) { t[0] = ComputeXIntercept(w, degree); return 1; } break; } } /* Otherwise, solve recursively after */ /* subdividing control polygon */ Bezier(w, degree, 0.5, Left, Right); left_count = FindRoots(Left, degree, left_t, depth+1); right_count = FindRoots(Right, degree, right_t, depth+1); /* Gather solutions together */ for (i = 0; i < left_count; i++) { t[i] = left_t[i]; } for (i = 0; i < right_count; i++) { t[i+left_count] = right_t[i]; } /* Send back total number of solutions */ return (left_count+right_count); } /* * CrossingCount : * Count the number of times a Bezier control polygon * crosses the 0-axis. This number is >= the number of roots. * */ static int CrossingCount(V, degree) Point2 *V; /* Control pts of Bezier curve */ int degree; /* Degreee of Bezier curve */ { int i; int n_crossings = 0; /* Number of zero-crossings */ int sign, old_sign; /* Sign of coefficients */ sign = old_sign = SGN(V[0].y); for (i = 1; i <= degree; i++) { sign = SGN(V[i].y); if (sign != old_sign) n_crossings++; old_sign = sign; } return n_crossings; } /* * ControlPolygonFlatEnough : * Check if the control polygon of a Bezier curve is flat enough * for recursive subdivision to bottom out. * */ static int ControlPolygonFlatEnough(V, degree) Point2 *V; /* Control points */ int degree; /* Degree of polynomial */ { int i; /* Index variable */ double *distance; /* Distances from pts to line */ double max_distance_above; /* maximum of these */ double max_distance_below; double error; /* Precision of root */ double intercept_1, intercept_2, left_intercept, right_intercept; double a, b, c; /* Coefficients of implicit */ /* eqn for line from V[0]-V[deg]*/ /* Find the perpendicular distance */ /* from each interior control point to */ /* line connecting V[0] and V[degree] */ distance = (double *)malloc((unsigned)(degree + 1) * sizeof(double)); { double abSquared; /* Derive the implicit equation for line connecting first *' /* and last control points */ a = V[0].y - V[degree].y; b = V[degree].x - V[0].x; c = V[0].x * V[degree].y - V[degree].x * V[0].y; abSquared = (a * a) + (b * b); for (i = 1; i < degree; i++) { /* Compute distance from each of the points to that line */ distance[i] = a * V[i].x + b * V[i].y + c; if (distance[i] > 0.0) { distance[i] = (distance[i] * distance[i]) / abSquared; } if (distance[i] < 0.0) { distance[i] = -((distance[i] * distance[i]) / abSquared); } } } /* Find the largest distance */ max_distance_above = 0.0; max_distance_below = 0.0; for (i = 1; i < degree; i++) { if (distance[i] < 0.0) { max_distance_below = MIN(max_distance_below, distance[i]); }; if (distance[i] > 0.0) { max_distance_above = MAX(max_distance_above, distance[i]); } } free((char *)distance); { double det, dInv; double a1, b1, c1, a2, b2, c2; /* Implicit equation for zero line */ a1 = 0.0; b1 = 1.0; c1 = 0.0; /* Implicit equation for "above" line */ a2 = a; b2 = b; c2 = c + max_distance_above; det = a1 * b2 - a2 * b1; dInv = 1.0/det; intercept_1 = (b1 * c2 - b2 * c1) * dInv; /* Implicit equation for "below" line */ a2 = a; b2 = b; c2 = c + max_distance_below; det = a1 * b2 - a2 * b1; dInv = 1.0/det; intercept_2 = (b1 * c2 - b2 * c1) * dInv; } /* Compute intercepts of bounding box */ left_intercept = MIN(intercept_1, intercept_2); right_intercept = MAX(intercept_1, intercept_2); error = 0.5 * (right_intercept-left_intercept); if (error < EPSILON) { return 1; } else { return 0; } } /* * ComputeXIntercept : * Compute intersection of chord from first control point to last * with 0-axis. * */ /* NOTE: "T" and "Y" do not have to be computed, and there are many useless * operations in the following (e.g. "0.0 - 0.0"). */ static double ComputeXIntercept(V, degree) Point2 *V; /* Control points */ int degree; /* Degree of curve */ { double XLK, YLK, XNM, YNM, XMK, YMK; double det, detInv; double S, T; double X, Y; XLK = 1.0 - 0.0; YLK = 0.0 - 0.0; XNM = V[degree].x - V[0].x; YNM = V[degree].y - V[0].y; XMK = V[0].x - 0.0; YMK = V[0].y - 0.0; det = XNM*YLK - YNM*XLK; detInv = 1.0/det; S = (XNM*YMK - YNM*XMK) * detInv; /* T = (XLK*YMK - YLK*XMK) * detInv; */ X = 0.0 + XLK * S; /* Y = 0.0 + YLK * S; */ return X; } /* * Bezier : * Evaluate a Bezier curve at a particular parameter value * Fill in control points for resulting sub-curves if "Left" and * "Right" are non-null. * */ static Point2 Bezier(V, degree, t, Left, Right) int degree; /* Degree of bezier curve */ Point2 *V; /* Control pts */ double t; /* Parameter value */ Point2 *Left; /* RETURN left half ctl pts */ Point2 *Right; /* RETURN right half ctl pts */ { int i, j; /* Index variables */ Point2 Vtemp[W_DEGREE+1][W_DEGREE+1]; /* Copy control points */ for (j =0; j <= degree; j++) { Vtemp[0][j] = V[j]; } /* Triangle computation */ for (i = 1; i <= degree; i++) { for (j =0 ; j <= degree - i; j++) { Vtemp[i][j].x = (1.0 - t) * Vtemp[i-1][j].x + t * Vtemp[i-1][j+1].x; Vtemp[i][j].y = (1.0 - t) * Vtemp[i-1][j].y + t * Vtemp[i-1][j+1].y; } } if (Left != NULL) { for (j = 0; j <= degree; j++) { Left[j] = Vtemp[j][0]; } } if (Right != NULL) { for (j = 0; j <= degree; j++) { Right[j] = Vtemp[degree-j][j]; } } return (Vtemp[degree][0]); } static Vector2 V2ScaleII(v, s) Vector2 *v; double s; { Vector2 result; result.x = v->x * s; result.y = v->y * s; return (result); }