src/third_party/skia/experimental/Intersection/CubicUtilities.cpp

   1 /*
   2  * Copyright 2012 Google Inc.
   3  *
   4  * Use of this source code is governed by a BSD-style license that can be
   5  * found in the LICENSE file.
   6  */
   7 #include "CubicUtilities.h"
   8 #include "Extrema.h"
   9 #include "LineUtilities.h"
  10 #include "QuadraticUtilities.h"
  11
  12 const int gPrecisionUnit = 256; // FIXME: arbitrary -- should try different values in test framework
  13
  14 // FIXME: cache keep the bounds and/or precision with the caller?
  15 double calcPrecision(const Cubic& cubic) {
  16     _Rect dRect;
  17     dRect.setBounds(cubic); // OPTIMIZATION: just use setRawBounds ?
  18     double width = dRect.right - dRect.left;
  19     double height = dRect.bottom - dRect.top;
  20     return (width > height ? width : height) / gPrecisionUnit;
  21 }
  22
  23 #ifdef SK_DEBUG
  24 double calcPrecision(const Cubic& cubic, double t, double scale) {
  25     Cubic part;
  26     sub_divide(cubic, SkTMax(0., t - scale), SkTMin(1., t + scale), part);
  27     return calcPrecision(part);
  28 }
  29 #endif
  30
  31 bool clockwise(const Cubic& c) {
  32     double sum = (c[0].x - c[3].x) * (c[0].y + c[3].y);
  33     for (int idx = 0; idx < 3; ++idx){
  34         sum += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y);
  35     }
  36     return sum <= 0;
  37 }
  38
  39 void coefficients(const double* cubic, double& A, double& B, double& C, double& D) {
  40     A = cubic[6]; // d
  41     B = cubic[4] * 3; // 3*c
  42     C = cubic[2] * 3; // 3*b
  43     D = cubic[0]; // a
  44     A -= D - C + B;     // A =   -a + 3*b - 3*c + d
  45     B += 3 * D - 2 * C; // B =  3*a - 6*b + 3*c
  46     C -= 3 * D;         // C = -3*a + 3*b
  47 }
  48
  49 bool controls_contained_by_ends(const Cubic& c) {
  50     _Vector startTan = c[1] - c[0];
  51     if (startTan.x == 0 && startTan.y == 0) {
  52         startTan = c[2] - c[0];
  53     }
  54     _Vector endTan = c[2] - c[3];
  55     if (endTan.x == 0 && endTan.y == 0) {
  56         endTan = c[1] - c[3];
  57     }
  58     if (startTan.dot(endTan) >= 0) {
  59         return false;
  60     }
  61     _Line startEdge = {c[0], c[0]};
  62     startEdge[1].x -= startTan.y;
  63     startEdge[1].y += startTan.x;
  64     _Line endEdge = {c[3], c[3]};
  65     endEdge[1].x -= endTan.y;
  66     endEdge[1].y += endTan.x;
  67     double leftStart1 = is_left(startEdge, c[1]);
  68     if (leftStart1 * is_left(startEdge, c[2]) < 0) {
  69         return false;
  70     }
  71     double leftEnd1 = is_left(endEdge, c[1]);
  72     if (leftEnd1 * is_left(endEdge, c[2]) < 0) {
  73         return false;
  74     }
  75     return leftStart1 * leftEnd1 >= 0;
  76 }
  77
  78 bool ends_are_extrema_in_x_or_y(const Cubic& c) {
  79     return (between(c[0].x, c[1].x, c[3].x) && between(c[0].x, c[2].x, c[3].x))
  80             || (between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].y));
  81 }
  82
  83 bool monotonic_in_y(const Cubic& c) {
  84     return between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].y);
  85 }
  86
  87 bool serpentine(const Cubic& c) {
  88     if (!controls_contained_by_ends(c)) {
  89         return false;
  90     }
  91     double wiggle = (c[0].x - c[2].x) * (c[0].y + c[2].y);
  92     for (int idx = 0; idx < 2; ++idx){
  93         wiggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y);
  94     }
  95     double waggle = (c[1].x - c[3].x) * (c[1].y + c[3].y);
  96     for (int idx = 1; idx < 3; ++idx){
  97         waggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y);
  98     }
  99     return wiggle * waggle < 0;
 100 }
 101
 102 // cubic roots
 103
 104 const double PI = 4 * atan(1);
 105
 106 // from SkGeometry.cpp (and Numeric Solutions, 5.6)
 107 int cubicRootsValidT(double A, double B, double C, double D, double t[3]) {
 108 #if 0
 109     if (approximately_zero(A)) {  // we're just a quadratic
 110         return quadraticRootsValidT(B, C, D, t);
 111     }
 112     double a, b, c;
 113     {
 114         double invA = 1 / A;
 115         a = B * invA;
 116         b = C * invA;
 117         c = D * invA;
 118     }
 119     double a2 = a * a;
 120     double Q = (a2 - b * 3) / 9;
 121     double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
 122     double Q3 = Q * Q * Q;
 123     double R2MinusQ3 = R * R - Q3;
 124     double adiv3 = a / 3;
 125     double* roots = t;
 126     double r;
 127
 128     if (R2MinusQ3 < 0)   // we have 3 real roots
 129     {
 130         double theta = acos(R / sqrt(Q3));
 131         double neg2RootQ = -2 * sqrt(Q);
 132
 133         r = neg2RootQ * cos(theta / 3) - adiv3;
 134         if (is_unit_interval(r))
 135             *roots++ = r;
 136
 137         r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
 138         if (is_unit_interval(r))
 139             *roots++ = r;
 140
 141         r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
 142         if (is_unit_interval(r))
 143             *roots++ = r;
 144     }
 145     else                // we have 1 real root
 146     {
 147         double A = fabs(R) + sqrt(R2MinusQ3);
 148         A = cube_root(A);
 149         if (R > 0) {
 150             A = -A;
 151         }
 152         if (A != 0) {
 153             A += Q / A;
 154         }
 155         r = A - adiv3;
 156         if (is_unit_interval(r))
 157             *roots++ = r;
 158     }
 159     return (int)(roots - t);
 160 #else
 161     double s[3];
 162     int realRoots = cubicRootsReal(A, B, C, D, s);
 163     int foundRoots = add_valid_ts(s, realRoots, t);
 164     return foundRoots;
 165 #endif
 166 }
 167
 168 int cubicRootsReal(double A, double B, double C, double D, double s[3]) {
 169 #ifdef SK_DEBUG
 170     // create a string mathematica understands
 171     // GDB set print repe 15 # if repeated digits is a bother
 172     //     set print elements 400 # if line doesn't fit
 173     char str[1024];
 174     bzero(str, sizeof(str));
 175     sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D);
 176     mathematica_ize(str, sizeof(str));
 177 #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA
 178     SkDebugf("%s\n", str);
 179 #endif
 180 #endif
 181     if (approximately_zero(A)
 182             && approximately_zero_when_compared_to(A, B)
 183             && approximately_zero_when_compared_to(A, C)
 184             && approximately_zero_when_compared_to(A, D)) {  // we're just a quadratic
 185         return quadraticRootsReal(B, C, D, s);
 186     }
 187     if (approximately_zero_when_compared_to(D, A)
 188             && approximately_zero_when_compared_to(D, B)
 189             && approximately_zero_when_compared_to(D, C)) { // 0 is one root
 190         int num = quadraticRootsReal(A, B, C, s);
 191         for (int i = 0; i < num; ++i) {
 192             if (approximately_zero(s[i])) {
 193                 return num;
 194             }
 195         }
 196         s[num++] = 0;
 197         return num;
 198     }
 199     if (approximately_zero(A + B + C + D)) { // 1 is one root
 200         int num = quadraticRootsReal(A, A + B, -D, s);
 201         for (int i = 0; i < num; ++i) {
 202             if (AlmostEqualUlps(s[i], 1)) {
 203                 return num;
 204             }
 205         }
 206         s[num++] = 1;
 207         return num;
 208     }
 209     double a, b, c;
 210     {
 211         double invA = 1 / A;
 212         a = B * invA;
 213         b = C * invA;
 214         c = D * invA;
 215     }
 216     double a2 = a * a;
 217     double Q = (a2 - b * 3) / 9;
 218     double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
 219     double R2 = R * R;
 220     double Q3 = Q * Q * Q;
 221     double R2MinusQ3 = R2 - Q3;
 222     double adiv3 = a / 3;
 223     double r;
 224     double* roots = s;
 225 #if 0
 226     if (approximately_zero_squared(R2MinusQ3) && AlmostEqualUlps(R2, Q3)) {
 227         if (approximately_zero_squared(R)) {/* one triple solution */
 228             *roots++ = -adiv3;
 229         } else { /* one single and one double solution */
 230
 231             double u = cube_root(-R);
 232             *roots++ = 2 * u - adiv3;
 233             *roots++ = -u - adiv3;
 234         }
 235     }
 236     else
 237 #endif
 238     if (R2MinusQ3 < 0)   // we have 3 real roots
 239     {
 240         double theta = acos(R / sqrt(Q3));
 241         double neg2RootQ = -2 * sqrt(Q);
 242
 243         r = neg2RootQ * cos(theta / 3) - adiv3;
 244         *roots++ = r;
 245
 246         r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
 247         if (!AlmostEqualUlps(s[0], r)) {
 248             *roots++ = r;
 249         }
 250         r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
 251         if (!AlmostEqualUlps(s[0], r) && (roots - s == 1 || !AlmostEqualUlps(s[1], r))) {
 252             *roots++ = r;
 253         }
 254     }
 255     else                // we have 1 real root
 256     {
 257         double sqrtR2MinusQ3 = sqrt(R2MinusQ3);
 258         double A = fabs(R) + sqrtR2MinusQ3;
 259         A = cube_root(A);
 260         if (R > 0) {
 261             A = -A;
 262         }
 263         if (A != 0) {
 264             A += Q / A;
 265         }
 266         r = A - adiv3;
 267         *roots++ = r;
 268         if (AlmostEqualUlps(R2, Q3)) {
 269             r = -A / 2 - adiv3;
 270             if (!AlmostEqualUlps(s[0], r)) {
 271                 *roots++ = r;
 272             }
 273         }
 274     }
 275     return (int)(roots - s);
 276 }
 277
 278 // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf
 279 // c(t)  = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3
 280 // c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2
 281 //       = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2
 282 static double derivativeAtT(const double* cubic, double t) {
 283     double one_t = 1 - t;
 284     double a = cubic[0];
 285     double b = cubic[2];
 286     double c = cubic[4];
 287     double d = cubic[6];
 288     return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t);
 289 }
 290
 291 double dx_at_t(const Cubic& cubic, double t) {
 292     return derivativeAtT(&cubic[0].x, t);
 293 }
 294
 295 double dy_at_t(const Cubic& cubic, double t) {
 296     return derivativeAtT(&cubic[0].y, t);
 297 }
 298
 299 // OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t?
 300 _Vector dxdy_at_t(const Cubic& cubic, double t) {
 301     _Vector result = { derivativeAtT(&cubic[0].x, t), derivativeAtT(&cubic[0].y, t) };
 302     return result;
 303 }
 304
 305 // OPTIMIZE? share code with formulate_F1DotF2
 306 int find_cubic_inflections(const Cubic& src, double tValues[])
 307 {
 308     double Ax = src[1].x - src[0].x;
 309     double Ay = src[1].y - src[0].y;
 310     double Bx = src[2].x - 2 * src[1].x + src[0].x;
 311     double By = src[2].y - 2 * src[1].y + src[0].y;
 312     double Cx = src[3].x + 3 * (src[1].x - src[2].x) - src[0].x;
 313     double Cy = src[3].y + 3 * (src[1].y - src[2].y) - src[0].y;
 314     return quadraticRootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues);
 315 }
 316
 317 static void formulate_F1DotF2(const double src[], double coeff[4])
 318 {
 319     double a = src[2] - src[0];
 320     double b = src[4] - 2 * src[2] + src[0];
 321     double c = src[6] + 3 * (src[2] - src[4]) - src[0];
 322     coeff[0] = c * c;
 323     coeff[1] = 3 * b * c;
 324     coeff[2] = 2 * b * b + c * a;
 325     coeff[3] = a * b;
 326 }
 327
 328 /*  from SkGeometry.cpp
 329     Looking for F' dot F'' == 0
 330
 331     A = b - a
 332     B = c - 2b + a
 333     C = d - 3c + 3b - a
 334
 335     F' = 3Ct^2 + 6Bt + 3A
 336     F'' = 6Ct + 6B
 337
 338     F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
 339 */
 340 int find_cubic_max_curvature(const Cubic& src, double tValues[])
 341 {
 342     double coeffX[4], coeffY[4];
 343     int i;
 344     formulate_F1DotF2(&src[0].x, coeffX);
 345     formulate_F1DotF2(&src[0].y, coeffY);
 346     for (i = 0; i < 4; i++) {
 347         coeffX[i] = coeffX[i] + coeffY[i];
 348     }
 349     return cubicRootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues);
 350 }
 351
 352
 353 bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) {
 354     double dy = cubic[index].y - cubic[zero].y;
 355     double dx = cubic[index].x - cubic[zero].x;
 356     if (approximately_zero(dy)) {
 357         if (approximately_zero(dx)) {
 358             return false;
 359         }
 360         memcpy(rotPath, cubic, sizeof(Cubic));
 361         return true;
 362     }
 363     for (int index = 0; index < 4; ++index) {
 364         rotPath[index].x = cubic[index].x * dx + cubic[index].y * dy;
 365         rotPath[index].y = cubic[index].y * dx - cubic[index].x * dy;
 366     }
 367     return true;
 368 }
 369
 370 #if 0 // unused for now
 371 double secondDerivativeAtT(const double* cubic, double t) {
 372     double a = cubic[0];
 373     double b = cubic[2];
 374     double c = cubic[4];
 375     double d = cubic[6];
 376     return (c - 2 * b + a) * (1 - t) + (d - 2 * c + b) * t;
 377 }
 378 #endif
 379
 380 _Point top(const Cubic& cubic, double startT, double endT) {
 381     Cubic sub;
 382     sub_divide(cubic, startT, endT, sub);
 383     _Point topPt = sub[0];
 384     if (topPt.y > sub[3].y || (topPt.y == sub[3].y && topPt.x > sub[3].x)) {
 385         topPt = sub[3];
 386     }
 387     double extremeTs[2];
 388     if (!monotonic_in_y(sub)) {
 389         int roots = findExtrema(sub[0].y, sub[1].y, sub[2].y, sub[3].y, extremeTs);
 390         for (int index = 0; index < roots; ++index) {
 391             _Point mid;
 392             double t = startT + (endT - startT) * extremeTs[index];
 393             xy_at_t(cubic, t, mid.x, mid.y);
 394             if (topPt.y > mid.y || (topPt.y == mid.y && topPt.x > mid.x)) {
 395                 topPt = mid;
 396             }
 397         }
 398     }
 399     return topPt;
 400 }
 401
 402 // OPTIMIZE: avoid computing the unused half
 403 void xy_at_t(const Cubic& cubic, double t, double& x, double& y) {
 404     _Point xy = xy_at_t(cubic, t);
 405     if (&x) {
 406         x = xy.x;
 407     }
 408     if (&y) {
 409         y = xy.y;
 410     }
 411 }
 412
 413 _Point xy_at_t(const Cubic& cubic, double t) {
 414     double one_t = 1 - t;
 415     double one_t2 = one_t * one_t;
 416     double a = one_t2 * one_t;
 417     double b = 3 * one_t2 * t;
 418     double t2 = t * t;
 419     double c = 3 * one_t * t2;
 420     double d = t2 * t;
 421     _Point result = {a * cubic[0].x + b * cubic[1].x + c * cubic[2].x + d * cubic[3].x,
 422             a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y};
 423     return result;
 424 }