1 /* 2 Solving the Nearest Point-on-Curve Problem 3 and 4 A Bezier Curve-Based Root-Finder 5 by Philip J. Schneider 6 from "Graphics Gems", Academic Press, 1990 7 */ 8 9 /* point_on_curve.c */ 10 11 #include <stdio.h> 12 #include <malloc.h> 13 #include <math.h> 14 #include "GraphicsGems.h" 15 16 #define TESTMODE 17 18 /* 19 * Forward declarations 20 */ 21 Point2 NearestPointOnCurve(); 22 static int FindRoots(); 23 static Point2 *ConvertToBezierForm(); 24 static double ComputeXIntercept(); 25 static int ControlPolygonFlatEnough(); 26 static int CrossingCount(); 27 static Point2 Bezier(); 28 static Vector2 V2ScaleII(); 29 30 int MAXDEPTH = 64; /* Maximum depth for recursion */ 31 32 #define EPSILON (ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */ 33 #define DEGREE 3 /* Cubic Bezier curve */ 34 #define W_DEGREE 5 /* Degree of eqn to find roots of */ 35 36 #ifdef TESTMODE 37 /* 38 * main : 39 * Given a cubic Bezier curve (i.e., its control points), and some 40 * arbitrary point in the plane, find the point on the curve 41 * closest to that arbitrary point. 42 */ 43 main() 44 { 45 46 static Point2 bezCurve[4] = { /* A cubic Bezier curve */ 47 { 0.0, 0.0 }, 48 { 1.0, 2.0 }, 49 { 3.0, 3.0 }, 50 { 4.0, 2.0 }, 51 }; 52 static Point2 arbPoint = { 3.5, 2.0 }; /*Some arbitrary point*/ 53 Point2 pointOnCurve; /* Nearest point on the curve */ 54 55 /* Find the closest point */ 56 pointOnCurve = NearestPointOnCurve(arbPoint, bezCurve); 57 printf("pointOnCurve : (%4.4f, %4.4f)\n", pointOnCurve.x, 58 pointOnCurve.y); 59 } 60 #endif /* TESTMODE */ 61 62 63 /* 64 * NearestPointOnCurve : 65 * Compute the parameter value of the point on a Bezier 66 * curve segment closest to some arbtitrary, user-input point. 67 * Return the point on the curve at that parameter value. 68 * 69 */ 70 Point2 NearestPointOnCurve(P, V) 71 Point2 P; /* The user-supplied point */ 72 Point2 *V; /* Control points of cubic Bezier */ 73 { 74 Point2 *w; /* Ctl pts for 5th-degree eqn */ 75 double t_candidate[W_DEGREE]; /* Possible roots */ 76 int n_solutions; /* Number of roots found */ 77 double t; /* Parameter value of closest pt*/ 78 79 /* Convert problem to 5th-degree Bezier form */ 80 w = ConvertToBezierForm(P, V); 81 82 /* Find all possible roots of 5th-degree equation */ 83 n_solutions = FindRoots(w, W_DEGREE, t_candidate, 0); 84 free((char *)w); 85 86 /* Compare distances of P to all candidates, and to t=0, and t=1 */ 87 { 88 double dist, new_dist; 89 Point2 p; 90 Vector2 v; 91 int i; 92 93 94 /* Check distance to beginning of curve, where t = 0 */ 95 dist = V2SquaredLength(V2Sub(&P, &V[0], &v)); 96 t = 0.0; 97 98 /* Find distances for candidate points */ 99 for (i = 0; i < n_solutions; i++) { 100 p = Bezier(V, DEGREE, t_candidate[i], 101 (Point2 *)NULL, (Point2 *)NULL); 102 new_dist = V2SquaredLength(V2Sub(&P, &p, &v)); 103 if (new_dist < dist) { 104 dist = new_dist; 105 t = t_candidate[i]; 106 } 107 } 108 109 /* Finally, look at distance to end point, where t = 1.0 */ 110 new_dist = V2SquaredLength(V2Sub(&P, &V[DEGREE], &v)); 111 if (new_dist < dist) { 112 dist = new_dist; 113 t = 1.0; 114 } 115 } 116 117 /* Return the point on the curve at parameter value t */ 118 printf("t : %4.12f\n", t); 119 return (Bezier(V, DEGREE, t, (Point2 *)NULL, (Point2 *)NULL)); 120 } 121 122 123 /* 124 * ConvertToBezierForm : 125 * Given a point and a Bezier curve, generate a 5th-degree 126 * Bezier-format equation whose solution finds the point on the 127 * curve nearest the user-defined point. 128 */ 129 static Point2 *ConvertToBezierForm(P, V) 130 Point2 P; /* The point to find t for */ 131 Point2 *V; /* The control points */ 132 { 133 int i, j, k, m, n, ub, lb; 134 int row, column; /* Table indices */ 135 Vector2 c[DEGREE+1]; /* V(i)'s - P */ 136 Vector2 d[DEGREE]; /* V(i+1) - V(i) */ 137 Point2 *w; /* Ctl pts of 5th-degree curve */ 138 double cdTable[3][4]; /* Dot product of c, d */ 139 static double z[3][4] = { /* Precomputed "z" for cubics */ 140 {1.0, 0.6, 0.3, 0.1}, 141 {0.4, 0.6, 0.6, 0.4}, 142 {0.1, 0.3, 0.6, 1.0}, 143 }; 144 145 146 /*Determine the c's -- these are vectors created by subtracting*/ 147 /* point P from each of the control points */ 148 for (i = 0; i <= DEGREE; i++) { 149 V2Sub(&V[i], &P, &c[i]); 150 } 151 /* Determine the d's -- these are vectors created by subtracting*/ 152 /* each control point from the next */ 153 for (i = 0; i <= DEGREE - 1; i++) { 154 d[i] = V2ScaleII(V2Sub(&V[i+1], &V[i], &d[i]), 3.0); 155 } 156 157 /* Create the c,d table -- this is a table of dot products of the */ 158 /* c's and d's */ 159 for (row = 0; row <= DEGREE - 1; row++) { 160 for (column = 0; column <= DEGREE; column++) { 161 cdTable[row][column] = V2Dot(&d[row], &c[column]); 162 } 163 } 164 165 /* Now, apply the z's to the dot products, on the skew diagonal*/ 166 /* Also, set up the x-values, making these "points" */ 167 w = (Point2 *)malloc((unsigned)(W_DEGREE+1) * sizeof(Point2)); 168 for (i = 0; i <= W_DEGREE; i++) { 169 w[i].y = 0.0; 170 w[i].x = (double)(i) / W_DEGREE; 171 } 172 173 n = DEGREE; 174 m = DEGREE-1; 175 for (k = 0; k <= n + m; k++) { 176 lb = MAX(0, k - m); 177 ub = MIN(k, n); 178 for (i = lb; i <= ub; i++) { 179 j = k - i; 180 w[i+j].y += cdTable[j][i] * z[j][i]; 181 } 182 } 183 184 return (w); 185 } 186 187 188 /* 189 * FindRoots : 190 * Given a 5th-degree equation in Bernstein-Bezier form, find 191 * all of the roots in the interval [0, 1]. Return the number 192 * of roots found. 193 */ 194 static int FindRoots(w, degree, t, depth) 195 Point2 *w; /* The control points */ 196 int degree; /* The degree of the polynomial */ 197 double *t; /* RETURN candidate t-values */ 198 int depth; /* The depth of the recursion */ 199 { 200 int i; 201 Point2 Left[W_DEGREE+1], /* New left and right */ 202 Right[W_DEGREE+1]; /* control polygons */ 203 int left_count, /* Solution count from */ 204 right_count; /* children */ 205 double left_t[W_DEGREE+1], /* Solutions from kids */ 206 right_t[W_DEGREE+1]; 207 208 switch (CrossingCount(w, degree)) { 209 case 0 : { /* No solutions here */ 210 return 0; 211 } 212 case 1 : { /* Unique solution */ 213 /* Stop recursion when the tree is deep enough */ 214 /* if deep enough, return 1 solution at midpoint */ 215 if (depth >= MAXDEPTH) { 216 t[0] = (w[0].x + w[W_DEGREE].x) / 2.0; 217 return 1; 218 } 219 if (ControlPolygonFlatEnough(w, degree)) { 220 t[0] = ComputeXIntercept(w, degree); 221 return 1; 222 } 223 break; 224 } 225 } 226 227 /* Otherwise, solve recursively after */ 228 /* subdividing control polygon */ 229 Bezier(w, degree, 0.5, Left, Right); 230 left_count = FindRoots(Left, degree, left_t, depth+1); 231 right_count = FindRoots(Right, degree, right_t, depth+1); 232 233 234 /* Gather solutions together */ 235 for (i = 0; i < left_count; i++) { 236 t[i] = left_t[i]; 237 } 238 for (i = 0; i < right_count; i++) { 239 t[i+left_count] = right_t[i]; 240 } 241 242 /* Send back total number of solutions */ 243 return (left_count+right_count); 244 } 245 246 247 /* 248 * CrossingCount : 249 * Count the number of times a Bezier control polygon 250 * crosses the 0-axis. This number is >= the number of roots. 251 * 252 */ 253 static int CrossingCount(V, degree) 254 Point2 *V; /* Control pts of Bezier curve */ 255 int degree; /* Degreee of Bezier curve */ 256 { 257 int i; 258 int n_crossings = 0; /* Number of zero-crossings */ 259 int sign, old_sign; /* Sign of coefficients */ 260 261 sign = old_sign = SGN(V[0].y); 262 for (i = 1; i <= degree; i++) { 263 sign = SGN(V[i].y); 264 if (sign != old_sign) n_crossings++; 265 old_sign = sign; 266 } 267 return n_crossings; 268 } 269 270 271 272 /* 273 * ControlPolygonFlatEnough : 274 * Check if the control polygon of a Bezier curve is flat enough 275 * for recursive subdivision to bottom out. 276 * 277 */ 278 static int ControlPolygonFlatEnough(V, degree) 279 Point2 *V; /* Control points */ 280 int degree; /* Degree of polynomial */ 281 { 282 int i; /* Index variable */ 283 double *distance; /* Distances from pts to line */ 284 double max_distance_above; /* maximum of these */ 285 double max_distance_below; 286 double error; /* Precision of root */ 287 double intercept_1, 288 intercept_2, 289 left_intercept, 290 right_intercept; 291 double a, b, c; /* Coefficients of implicit */ 292 /* eqn for line from V[0]-V[deg]*/ 293 294 /* Find the perpendicular distance */ 295 /* from each interior control point to */ 296 /* line connecting V[0] and V[degree] */ 297 distance = (double *)malloc((unsigned)(degree + 1) * sizeof(double)); 298 { 299 double abSquared; 300 301 /* Derive the implicit equation for line connecting first *' 302 /* and last control points */ 303 a = V[0].y - V[degree].y; 304 b = V[degree].x - V[0].x; 305 c = V[0].x * V[degree].y - V[degree].x * V[0].y; 306 307 abSquared = (a * a) + (b * b); 308 309 for (i = 1; i < degree; i++) { 310 /* Compute distance from each of the points to that line */ 311 distance[i] = a * V[i].x + b * V[i].y + c; 312 if (distance[i] > 0.0) { 313 distance[i] = (distance[i] * distance[i]) / abSquared; 314 } 315 if (distance[i] < 0.0) { 316 distance[i] = -((distance[i] * distance[i]) / abSquared); 317 } 318 } 319 } 320 321 322 /* Find the largest distance */ 323 max_distance_above = 0.0; 324 max_distance_below = 0.0; 325 for (i = 1; i < degree; i++) { 326 if (distance[i] < 0.0) { 327 max_distance_below = MIN(max_distance_below, distance[i]); 328 }; 329 if (distance[i] > 0.0) { 330 max_distance_above = MAX(max_distance_above, distance[i]); 331 } 332 } 333 free((char *)distance); 334 335 { 336 double det, dInv; 337 double a1, b1, c1, a2, b2, c2; 338 339 /* Implicit equation for zero line */ 340 a1 = 0.0; 341 b1 = 1.0; 342 c1 = 0.0; 343 344 /* Implicit equation for "above" line */ 345 a2 = a; 346 b2 = b; 347 c2 = c + max_distance_above; 348 349 det = a1 * b2 - a2 * b1; 350 dInv = 1.0/det; 351 352 intercept_1 = (b1 * c2 - b2 * c1) * dInv; 353 354 /* Implicit equation for "below" line */ 355 a2 = a; 356 b2 = b; 357 c2 = c + max_distance_below; 358 359 det = a1 * b2 - a2 * b1; 360 dInv = 1.0/det; 361 362 intercept_2 = (b1 * c2 - b2 * c1) * dInv; 363 } 364 365 /* Compute intercepts of bounding box */ 366 left_intercept = MIN(intercept_1, intercept_2); 367 right_intercept = MAX(intercept_1, intercept_2); 368 369 error = 0.5 * (right_intercept-left_intercept); 370 if (error < EPSILON) { 371 return 1; 372 } 373 else { 374 return 0; 375 } 376 } 377 378 379 380 /* 381 * ComputeXIntercept : 382 * Compute intersection of chord from first control point to last 383 * with 0-axis. 384 * 385 */ 386 /* NOTE: "T" and "Y" do not have to be computed, and there are many useless 387 * operations in the following (e.g. "0.0 - 0.0"). 388 */ 389 static double ComputeXIntercept(V, degree) 390 Point2 *V; /* Control points */ 391 int degree; /* Degree of curve */ 392 { 393 double XLK, YLK, XNM, YNM, XMK, YMK; 394 double det, detInv; 395 double S, T; 396 double X, Y; 397 398 XLK = 1.0 - 0.0; 399 YLK = 0.0 - 0.0; 400 XNM = V[degree].x - V[0].x; 401 YNM = V[degree].y - V[0].y; 402 XMK = V[0].x - 0.0; 403 YMK = V[0].y - 0.0; 404 405 det = XNM*YLK - YNM*XLK; 406 detInv = 1.0/det; 407 408 S = (XNM*YMK - YNM*XMK) * detInv; 409 /* T = (XLK*YMK - YLK*XMK) * detInv; */ 410 411 X = 0.0 + XLK * S; 412 /* Y = 0.0 + YLK * S; */ 413 414 return X; 415 } 416 417 418 /* 419 * Bezier : 420 * Evaluate a Bezier curve at a particular parameter value 421 * Fill in control points for resulting sub-curves if "Left" and 422 * "Right" are non-null. 423 * 424 */ 425 static Point2 Bezier(V, degree, t, Left, Right) 426 int degree; /* Degree of bezier curve */ 427 Point2 *V; /* Control pts */ 428 double t; /* Parameter value */ 429 Point2 *Left; /* RETURN left half ctl pts */ 430 Point2 *Right; /* RETURN right half ctl pts */ 431 { 432 int i, j; /* Index variables */ 433 Point2 Vtemp[W_DEGREE+1][W_DEGREE+1]; 434 435 436 /* Copy control points */ 437 for (j =0; j <= degree; j++) { 438 Vtemp[0][j] = V[j]; 439 } 440 441 /* Triangle computation */ 442 for (i = 1; i <= degree; i++) { 443 for (j =0 ; j <= degree - i; j++) { 444 Vtemp[i][j].x = 445 (1.0 - t) * Vtemp[i-1][j].x + t * Vtemp[i-1][j+1].x; 446 Vtemp[i][j].y = 447 (1.0 - t) * Vtemp[i-1][j].y + t * Vtemp[i-1][j+1].y; 448 } 449 } 450 451 if (Left != NULL) { 452 for (j = 0; j <= degree; j++) { 453 Left[j] = Vtemp[j][0]; 454 } 455 } 456 if (Right != NULL) { 457 for (j = 0; j <= degree; j++) { 458 Right[j] = Vtemp[degree-j][j]; 459 } 460 } 461 462 return (Vtemp[degree][0]); 463 } 464 465 static Vector2 V2ScaleII(v, s) 466 Vector2 *v; 467 double s; 468 { 469 Vector2 result; 470 471 result.x = v->x * s; result.y = v->y * s; 472 return (result); 473 } 474
