| OLD | NEW |
|---|
| (Empty) |
| 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) * siz
eof(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 } | |
| OLD | NEW |
|---|
Chromium Code Reviews
Issue 867213004:
remove prototype pathops code (Closed)
Base URL: https://skia.googlesource.com/skia.git@master